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ABSTRACT 

The  capacitated  generalized  transshipment  problem  is  the  most  general 
and  universally  applicable  member  of  the  class  of  network  optimization 
models.  This  model  subsumes,  as  specializations,  the  capacitated  and 
uncapacitated  transportation  problems  as  well  as  the  pure  network  special- 
izations of  these  models,  which  include  the  personnel  assignment  oroblem, 
the  maximum  flow,  and  shortest  path  formulations.  The  generalized 
network  problem,  in  turn,  can  be  viewed  as  a  specialization  of  a  linear 
programming  problem  having  at  most  two  non-zero  entries  in  each  column  of 
the  constraint  matrix.  A  detailed  description  is  given  of  the  implementa- 
tion of  an  efficient  algorithm  and  its  supporting  data  structures,  used 
to  solve  large-scale,  minimum-cost  generalized  transshipment  problems  on 
an  Apple  II  (64K)  microcomputer.  A  suite  of  advanced  techniques  for 
managing  minimum-cost  network  flow  models  and  inherent  data  elements  will 
also  be  discussed. 
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NOTATION 

Except  as  otherwise  indicated  in  the  text,  the  following  notational 
conventions  have  been  used  for  all  mathematical  expressions: 

Vectors:  Lower  case  Latin  (e.g.,  u) 

Vector  Components:     Lower  case  Latin  with  subscript  (e.g.,  u-) 

Matrices:  Upper  case  Latin  (e.g..  A) 

Matrix  Column:        Upper  case  Latin  with  superscript  column  index 

(e.g.,  N^ 
Matrix  Row:  Upper  case  Latin  with  subscript  row  index 

(e.g..  A.) 
Matrix  Element:       Lower  case  Latin  with  subscript  row  and  column 

indices  (e.g. ,  a. .) 
Set:  Upper  case  script  (e.g.,  ar) 

Scalars:  Lower  case  script  (e.g.,  g)    if  emphasis  is  required; 

otherwise  lower  case  Latin  (e.g.,  i) 


I.   INTRODUCTION 

Since  the  development  of  the  Simplex  method  by  George  Dantzig,  and 
the  introduction  of  the  transportation  model  by  Tjalling  Koopmans  (both 
circa  1947),  network  models  have  enjoyed  wide  use.  Perhaps  two  reasons 
for  this  attention  are  the  frequent  occurrence  of  situations  which  are 
readily  modelled  as  networks  and  the  mathematical  and  computational 
elegance  which  may  be  achieved  through  network  specializations  of  the 
Simplex  method.  Undoubtedly,  the  visually  appealing  graphical  description 
provided  by  the  network  formulation  has  contributed  much  to  the  managerial 
acceptance  of  these  models.  Network  models  have  been  used  in  a  large 
number  of  applications.  Jensen  and  Barnes  [Ref.  1]  provide  a  number  of 
examples  of  network  modelling  techniques  and  applications,  as  do  Kenninqton 
and  Helgason  [Ref.  2],  Dantzig  [Ref.  3],  and  Bradley  [Ref.  4].  Some  of 
these  applications  deal  with  military  logistic  and  distribution  systems, 
communications,  and  pipeline  systems,  personnel  and  resource  assignments, 
and  production  planning. 

In  recent  years,  advances  in  solid  state  technology  have  enabled 
design  and  production  of  extremely  powerful  microprocessor-based  computers. 
The  usefulness  of  these  computers  to  applications  of  mathematical  orogram- 
ming  has  been  largely  overlooked  by  all  but  a  few  researchers.  The 
computational  efficiency,  speed,  and  elegance  of  network  algorithms  and 
the  broad  range  of  application  of  the  generalized  network  formulation 
make  microcomputer-based  network  optimization  extremely  attractive. 


The  first  microcomputer-based  network  software  suite  has  been  designed 
and  implemented  on  an  Apple  II  Plus.  This  network  system  exhibits  more 
capability  than  any  other  package  available  (on  any  host  computer)  at 
this  writing.   It  is  designed  as  an  integrated  suite  of  programs  capable 
of  solving  pure  capacitated,  non-linear,  elastic  (with  fixed  charges), 
and  capacitated  generalized  network  problems  and  includes  a  user-friendly 
interface  which  facilitates  data  input  and  manipulation. 

The  capacitated  generalized  transshipment  problem  is  the  most  general 
and  universally  applicable  of  the  network  optimization  models.  This 
model  class  embraces,  as  specializations,  the  capacitated  transportation 
problem  and  the  pure  network  class  of  models.  As  viewed  here,  the  object 
of  these  formulations  is  to  determine  in  what  manner  a  good  or  commodity 
should  flow  through  a  network  such  that  flow  is  conserved  at  each  node 
and  the  total  cost  of  flow  through  the  network  is  minimized. 

(LP)    minimize  f(x)  cost 

s.t.     Ax  =  b         conservation  constraints 
lb  <_  x  <_  cp     bounds  on  flow 

This  problem  can  easily  be  formulated  as  a  linear  program  (LP),  however, 
the  network  specialization  provides  significant  computational  savings, 
often  producing  solutions  in  one  hundredth  the  time  [Ref.  2]  required  by 
the  equivalent  linear  program.  Additionally,  the  network  formulation, 
when  viewed  pictorially  as  a  collection  of  arcs  and  nodes,  has  an  obvious 
intuitive  appeal,  making  it  more  readily  accepted  and  understood  by 
non-analysts  [Ref.  4]. 
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There  is  an  important  difference  between  the  arcs  of  a  pure  network 
and  the  arcs  of  a  generalized  network;  associated  with  each  arc  of  a 
generalized  network  is  a  multiplier  which  acts  on  the  flow  through  that 
arc.  The  arc  multipliers  may  serve  to  transform  the  units  of  flow  or 
they  may  change  the  amount  of  flow  through  an  arc  [e.g.,  Ref.  5].  For 
example,  if  we  wish  to  represent  the  conversion  of  steel  into  automobile 
chasses  at  the  rate  of  1/10  ton  of  steel  per  chassis,  the  arc  multiplier 
converting  tons  of  steel  into  chasses  would  be  10.  Ten  automobile 
chasses  can  be  manufactured  from  each  ton  of  steel.  Alternatively,  if 
flow  on  an  arc  is  in  terms  of  investment  dollars  from  one  year  to  the 
next  at  an  annual  rate  of  return  of  12  percent,  the  multiplier  associated 
with  the  arc  representing  that  investment  would  be  1.12. 

In  the  spirit  of  Bradley,  Brown,  and  Graves  [Ref.  5]  and  Brown  and 
McBride  [Ref.  7],  the  network  formulation  may  be  described  as  a  directed 
graph  G  defined  by  a  set  of  nodes  ND   and  a  set  of  arcs  AR.     Henceforward, 
n  will  be  referred  to  as  the  number  of  arcs  in  a  network  and  m   will 
represent  the  number  of  nodes.  The  conceptual  constraint  matrix  A  is 
thus  m   rows  by  n   columns. 

Members  of  the  arc  set  are  indexed  by  k  and  are  defined  as  an  ordered 
pair  (tail,  head)  or  (source  node,  destination  node).  Associated  with 
each  arc  there  is  a  cost  per  unit  flow  c.  ,  a  lower  bound  on  flow  lb,  , 
an  upper  bound  on  flow,  or  capacity,  cp.  .  The  flow  on  arc  k  is 
represented  by  x,  . 

The  generalized  network  model  employs  an  arc  multiplier  m.  which 
represents  a  gain  or  loss  in  material  flowing  across  arc  k.  When  this 
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arc  multiplier  is  unity  for  all  arcs  in  the  network,  the  model  is  then  a 
pure  network  specialization  of  the  generalized  network  model. 

Each  node  can  be  designated  as  a  supply  node  (material  enters  the 
network),  a  demand  node  (material  leaves  the  network),  or  a  transshipment 
node  (material  just  passes  through). 

The  problem  is  to  minimize  the  total  cost  of  flow  on  all  the  arcs, 
such -that  the  flow  on  each  of  the  arcs  remains  within  the  stipulated 
bounds,  demands  are  met  from  available  supplies,  and  flow  is  conserved  at 
each  node.  Conservation  of  flow  requires  that  the  flow  leaving  a  node 
minus  the  flow  entering  a  node  equals  the  external  flow  or  requirement  of 
that  node.  In  generalized  networks,  the  flow  entering  a  node  is  the  sum 
of  the  flows  on  the  incoming  arcs  multiplied  by  the  associated  arc  gains. 
These  requirements  can  be  written  as: 

(GNP)    min    J^  c^x,^ 
ke^i? 

s.t.         2     X.    -         J2      \^](  ~  '^i»   ^"^    ^  ^       conservation 
keAff  'KtAR  of  flow 

leaving   i      entering   i 

lb  <_  X  <_  cp,  bounds  on  flow 

where  0  <  b.  =  supply  quantity,  if  i  is  a  supply  node 

0  >  b.  =  -  (demand  quantity),  if  i  is  a  demand  node 

0=b.  =0,  ifiisa  transshipment  node. 

The  convention  followed  here  is  that  supplies  are  represented  bv 
positive  magnitudes  and  demands  are  negative.  This  algebraic  template 
can  accommodate  a  model  with  inequality  flow  constraints,  insuring  that 
no  more  than  the  available  supply  will  be  utilized  and  that  no  less  than 
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the  demand  will  be  provided.  Slack  arcs  represent  the  difference  between 
available  supply  and  that  portion  actually  used,  and  surplus  arcs  measure 
the  extent  to  which  shipments  exceed  demand.  When  slack  and  surplus  arcs 
are  utilized  in  this  manner,  a  formulation  which  uses  inequality  con- 
straints can  be  transformed  to  the  equality  constraints  of  (GNP).  Slack 
and  surplus  arcs  are  considered  to  be  "logical"  arcs  as  opposed  to  the 
"structural"  arcs  of  the  original  problem. 

(GNP)  is  a  specialization  of  (LP).   If  each  column  of  the  constraint 
matrix  A  in  (LP)  corresponds  to  an  arc,  then  it  has  at  most  two  non-zero 
entries;  those  entries  can  be  scaled  to  be  +1  (for  the  tail)  and  some 
other  number  -m.  (at  the  head).  Thus,  the  constraint  matrix  of  (GNP) 
contains  elements  that  are  either  0,  1  or  -m,  (each  m.  is  admissably 
distinct).  When  the  Ax  matrix  multiplication  of  (LP)  is  enforced,  we 
obtain  the  flow  conservation  constraint  found  in  (GNP).  There  is  thus 
one  conservation  of  flow  constraint  in  (LP)  for  every  node  i:  A.x, 
where  A.  is  the  row  in  A  corresponding  to  the  i   node  in  (GNP), 
having  row  entries  of  +1  for  every  arc  originating  at  node  i,  and  -m. 
for  every  arc  terminating  at  node  i.  There  is  also  a  pair  of  bounds  in 
(LP)  and  in  (GNP)  for  eyery   arc  in  the  network. 

(GNP)  is  sufficiently  broad  to  enable  any  linear  program,  with  at 
most  two  non-zero  coefficients  associated  with  each  variable  (column),  to 
be  treated  as  a  generalized  network.  Even  when  the  linear  program  is  not 
entirely  composed  of  network  columns.  Brown  and  Wright  [Ref.  3]  have 
shown  that  many  real-life  linear  programs  contain  a  large  embedded 
network  structure  which,  once  discovered,  can  be  exploited  to  improve 
solution  efficiency  [Ref.  7]. 
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Solution  speed  and  efficient  data  storage  techniques  are  two  of  the 
most  attractive  features  of  the  network  model.  Due  to  the  obvious 
sparseness  of  the  network  constraint  matrix,  the  use  of  a  node-arc 
incidence  matrix  is  not  a  viable  method  of  basis  representation  for  data 
manipulation.  Space  and  speed-efficient  methods  of  handling  the  represen- 
tation and  manipulation  of  the  generalized  network  problem  are  introduced 
in  the  following  sections.  The  algorithm  which  will  be  described  is 
called  GENNET  [Ref.  7]. 
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II.  BASIS  REPRESENTATION  AND  DATA  STRUCTURES 

The  constraint  matrix  A  of  (LP)  represents  a  set  of  m   linear  equations 
in  n  unknowns.  If  m  <  n,   a  feasible  solution  to  these  constraints  may  be 
found  by  identifying  m   linearly  independent  columns  of  the  constraint 
matrix  A.  If  the  variables  associated  with  the  remaining  n-m  columns  of 
A  are  set  to  zero,  the  values  of  the  variables  corresponding  to  the  m 
linearly  independent  columns  may  be  uniquely  found  by  solving  the  result- 
ing set  of  exactly  determined  simultaneous  equations.  Any  set  of  m 
linearly  independent  columns  of  A  is  referred  to  as  a  basis.  The  solution 
to  the  set  of  linear  equations  obtained  by  setting  to  zero  the  variables 
corresponding  to  the  n-m   columns  of  A  not  included  in  the  basis  (i.e., 
non-basic  columns),  is  called  a  basic  solution. 

A  unique  characterizing  feature  of  the  pure  transshipment  problem  is 
the  triangular  nature  of  its  bases  [e.g.,  Ref.  3].  Triangular  bases  are 
particularly  convenient  because  they  are  easily  solved  by  direct  substitu- 
tion [e.g.,  Ref.  9].  The  solution  of  the  generalized  transshipment 
problem  is  somewhat  more  difficult  because  of  a  complication  introduced 
by  flow  mult i pi iers . 

What  is  now  known  as  the  capacitated  transshipment  problem  was  first 
posed,  with  a  discussion  of  solution  methods,  by  Koopmans  [Ref.  10]. 
Oantzig  [Ref.  3]  provides  the  first  general  exposition  of  solution 
technique  and  basis  representation  of  the  capacitated  generalized  trans- 
shipment problem,  referring  to  it  as  the  "Weighted  Distribution  Problem." 
Dantzig  shows  that  the  basis  of  a  generalized  transshipment  problem  has 
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unique  qualities  similar  to  the  basis  of  the  pure  network.  These  qual- 
ities will  be  exploited  in  the  problem  representation  and  solution 
approach  developed  here. 

A  graphical  representation  of  a  generalized  network  basis  results  in 
a  familiar  node  and  arc  display  [e.g.,  Ref.  2].  The  graph  thus  obtained 
is  not  necessarily  a  tree,  as  in  the  case  of  pure  networks,  but  may  be  a 
forest  whose  members  are  either  trees  or  trees  with  one  cycle  (or  loop). 
If  slack  variables  are  admissible,  then  the  network  must  also  include 
slack  arcs,  each  incident  with  only  one  node.  Figure  1  displays  a  node 
with  a  slack  arc. 


n 


Figure  1.  Node  with  a  Slack  Arc 

Danzig  [Ref.  3]  shows  that  each  component  of  the  basis  is  either  a 
rooted  tree  or  a  subgraph  with  one  cycle.  Figure  2  displays  typical 
components  of  a  generalized  network  basis.  A  rooted  tree  nay  be  v'ewed 
as  a  suograph  with  a  slack  arc  (and  thus  one  cycle)  at  the  root.  Each 
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Figure  2.  Typical  Components  of  a  Generalized  Network  3as1s 
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component  of  the  basis  is  a  connected  subgraph  with  precisely  one  cycle. 
The  basis  representation  of  a  generalized  network  is  "block  diagonal,"  as 
pictured  below.  The  diagonal  entries  are  submatrices  corresponding  to 
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the  component  trees/subgraphs  of  the  generalized  network  basis.  Those 
elements  representing  rooted  trees  are  upper  triangular,  as  they  are   in 
the  pure  network  case.  Oantzig  [Ref.  3]  shows  that  basis  components 
corresponding  to  subgraphs  having  one  cycle  may  be  put  in  "nearly  tri- 
angular" form,  with  only  one  column  having  a  non-zero  term  remaining 
below  the  diagonal.   If  the  variable  associated  with  that  "peculiar" 
column  is  treated  as  a  parameter,  values  for  the  variables  associated 
with  the  other  columns  can  be  obtained  in  terms  of  that  parameter.  These 
expressions  (in  terms  of  the  one  unknown  variable)  can  be  uniquely 
solved,  thereby  obtaining  a  complete  basic  solution  [Ref.  3].  A  solution 
method  for  generalized  networks  suggests  itself;  exploit  the  triangularity 
of  rooted  basis  components  with  efficient  data  structures  and  use  the 
method  proposed  by  Dantzig  for  nearly  triangular  basis  elements. 
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To  illustrate  these  concepts  suppose  we  have  the  basis  representation 
displ ayed  in  Figure  3. 


rooted  tree  with  slocK  arc  at  root 

Subgraph  with  one  cycle 

Figure  3.  Generalized  Network  Basis  Representation 

There  are  six  nodes  and  thus  six  arcs  in  this  basis,  one  of  which  is 
a  slack  arc  (arc  a)  corresponding  to  the  root  node  of  its  basis  tree. 
The  matrix  representation  of  this  basis  is  shown  in  Figure  4. 

Arc     a    b    c    d    e    f 
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Figure  4.   Matrix  Representation  of  a  Basis 
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The  multiplier  associated  with  slack  and  artificial  arcs  is  -1  and 

a  multiplier  of  +1  is'  used  for  surplus  arcs.  Thus,  m,  =  -1  since  arc 

a 

"a"  is  a  slack  arc.  Permuting  the  rows  of  this  basis  matrix  produces: 
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Figure  5.  Triangular/Nearly  Triangular  Basis 

The  dotted  lines  segregate  the  two  basis  components  B  and  B  . 
Note  that  B  corresponds  to  a  rooted  tree  and  is  thus  upper  triangular. 
B  is  "nearly  triangular."  There  is  one  "singleton"  column  in  the 
basis  which  corresponds  to  the  slack  arc  a. 

The  sparsity  of  the  constraint  matrix  and  the  graphical  structure 
of  the  model  lead  us  to  adopt  one-dimensional  data  structures  to  represent 
the  problem.  These  data  structures  provide  economical  storage  and  reduce 
computational  overhead.  The  two  most  common  types  of  data  structures 
found  in  network  optimization  algorithms  are  the  triple-lable  scheme, 
first  proposed  by  Ellis  Johnson  [Ref.  11]  and  the  preorder  traversal 
method,  as  described  by  Bradley,  Brown,  and  Graves  [Ref.  61. 
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The  triple-label  scheme  employs  three  functions  to  represent  the 
basis  graph:  a  predecessor  function  which  provides  the  parent,  father, 
or  immediate  ancestor  of  a  node;  a  successor  function  which  provides  the 
left-most  child  of  a  node  (sometimes  called  the  eldest  son);  and  a 
brother  function  which  provides  the  next  left-most  node  with  the  same 
parent  (i.e.,  the  next  oldest  brother  or  sibling).  These  functions  are 
exhibited  in  Figure  6. 
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Figure  5.   Triole-Label  Scheme 

Tne  triple-label  scheme  has  been  adootad  by  several  researcners 
[Ref.  12  for  the  (pure)  transportation  network  problem  and  Ref .    13  "*or 
generalized  networks].  Brown  and  McSride  [Ref.  7]  have  tested,  but 
not  adopted,  this  data  structure.  Kennington  and  Helgason  [Ref.  2] 
and  Jensen  and  Barnes  [Ref.    1]  repeat  textbook  explanations  of  the 
required  basis  update  actions  for  the  triple-label  iiethod  of  basis 
representation . 
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The  preorder  traversal  method  uses  two  functions  to  represent  and 
update  the  basis:  a  predecessor  function  P  and  a  preorder  function  IT. 
The  predecessor  function  provides  the  same  information  as  in  the  triple- 
label  scheme,  i.e.,  the  "father"  node.  The  predecessor  function  can  be 
easily  constructed  from  the  matrix  representation  of  the  basis  (e.g.. 
Figure  5).  After  rearranging  the  rows  and  columns  and  forming  triangular/ 
nearly  triangular  components  of  the  basis,  to  find  the  predecessor  of 
node  i : 

1.  Determine  the  row  number  of  node  i,  call  that  row  rl. 

2.  Determine  the  row  of  the  non-zero  off-diagonal  element  of  column 
rl,  cal 1  that  row  r2. 

3.  The  predecessor  of  node  i  is  the  node  represented  by  row  r2 .  If 

no  off-diagonal  element  is  found  in  step  2,  the  predecessor  of  node 
i  is  null.  (An  array  P  may  be  used  to  represent  the  predecessor 
function  and  a  distinguished  value,  e.g.,  m  +  I,   can  indicate  nodes 
with  no  predecessors.) 

Each  node  and  its  predecessor  define  the  basic  arc  (ignoring  for  the 
moment  the  arc's  orientation).  The  predecessor  function  can  be  used  to 
construct  a  "backpath"  from  any  node  to  the  root/cycle  of  its  basis 
component.  The  orientation  of  the  arc  in  the  original  network  can  be 
recorded  using  the  sign  of  the  predecessor.  If  the  arc  is  oriented  from 
node  i  to  P(i),  then  P(i)  >  0;  if  P(i)  <  0,  then  the  arc  is  directed  from 
P(i)  to  i. 

When  determining  the  predecessor  of  a  node  i,  if  a  multiplier  (-m|^) 
is  discovered  on  the  diagonal  of  the  nearly  triangular  basis  matrix,  then 
P(i)  <  0.  The  multiplier  value  is  associated  with  the  "destination  end" 
of  a  basic  arc.  The  arc  k  represented  by  a  column  having  -m^   on  the 
diagonal  of  the  basis  matrix  is  oriented  from  P(i)  to  i. 
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The  sign  of  the  predecessor  uniquely  determines  the  diagonal  entries 
of  the  basis  matrix.  If  P(i)  <  0,  the  diagonal  entry  in  the  row  represent- 
ing node  i  is  -m.  ;  if  P(i)  >  0,  the  diagonal  entry  is  +1.  A  predecessor 
function  applied  to  Figure  3  yields: 
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Figure  7.  Predecessor  Function 


A  node  may  have  "successors"  or  nodes  which  are  further  from  the 
root/cycle  than  the  node  in  question.  The  set  of  nodes  which  are  first 
encountered  on  all  paths  from  a  node,  except  the  path  to  the  root/cycle, 
are  called  the  "immediate  successors"  of  the  node.  A  basis  may  alterna- 
tively be  completely  represented  by  this  successor  rel ationshio.  However, 
because  any  node  can  have  a  variable  number  of  successors  but  only  one 
predecessor  in  the  basis,  the  predecessor  function  provides  a  more 
tractable  data  structure. 

In  the  example  presented  in  Figures  5  and  7,  node  4  is  the  predecessor 
of  both  nodes  5  and  6.  Clearly,  the  predecessor  function  does  not 
completely  represent  the  triangulated  basis.  The  piece  of  missing 
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information  is  whether  node  5  or  node  6  is  encountered  first  (top  to 
bottom)  in  the  triangulated  basis. 

The  technique  employed  here  to  represent  this  ordering  is  an  extension 
to  m-trees  of  the  preorder  binary  tree  traversal  described  by  Knuth  [Ref.  14] 
Graphically,  preorder  is  a  dynastic  ordering  reminiscent  of  the  inheritance 
of  thrones,  in  which  the  root  or  top-most  node  is  listed  first  and  its 
subtrees  are  listed  in  preorder,  until  every  node  is  listed.  Figure  3a 
displays  such  a  preorder. 
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Figure  3a.  Preorder  of  a  Tree 

Figure  3b  illustrates  an  extension  of  preorder  to  the  bases  of  general- 
ized networ*<3.  The  extension  of  oreorder  to  "trees"  with  cyclic  '^coxs  renuires 
that  one  noae  on  the  cycle  be  distinguished  as  the  first  oreorder  node.  The 
choice  of  that  distinguished  node  is  arbitrary  amongst  the  cycle  nodes. 
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Figure  3b.  extension  of  Preorder  to  Generalized  Networks 
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The  matrix  structure  of  the  triangulated  basis  also  induces  a  preorder. 
Row  i  always  precedes  its  successors  in  the  triangulated/nearly  tri- 
angulated basis  and  all  of  its  successors  precede  any  row  which  does  not 
precede  row  i  [Ref.  15]. 

To  summarize,  the  basis  representation  we  have  established  is  contained 
in  the  functions  P  and  IT.  P  indicates  the  predecessor  of  e^'jery   node  in 
the  generalized  network.  IT  provides  the  preorder-successor  of  each 
node.  By  "iterating"  IT  (e.g.,  IT(i),  IT(IT(i)),  IT( IT( IT( i ) ) )  ...  ), 
all  the  preorder-successors  of  node  i  may  eventually  be  found.  Iteration 
of  the  predecessor  and  preorder  functions  provides  a  means  of  tree 
traversal.  The  preorder  traversal  scheme  is  considered  to  be  the  more 
elegant  and  efficient  means  of  basis  representation  and  as  such  will  be 
used  in  this  description  of  the  GENNET  algorithm. 

The  predecessor  and  preorder  functions  are  implemented  as  arrays  P 
and  IT.  P(i)  indicates  the  predecessor  of  node  i.  Similarly,  IT(j) 
indicates  the  preorder-successor  of  node  j.  Associated  with  each  node  j 
is  an  arc  which  connects  node  j  with  its  predecessor  P(j).  This  arc 
corresponds  to  a  basic  variable.  The  index  of,  or  pointer  to,  the  arc 
connecting  node  j  and  node  P(j)  is  maintained  in  IVAR(j),  an  element  of  a 
node-length  (i.e.,  m   x  1)  array  IVAR.  The  sign  bit  of  P(j)  records  the 
orientation  of  the  arc  connecting  nodes  j  and  ?(j).   If  P(j)  >  0,  then 
the  arc  connecting  node  j  and  its  predecessor  is  oriented  from  node  ,i  to 
P(j);  an  arc  oriented  from  the  predecessor  of  node  j  to  node  j  is 
indicated  by  P(j)  <  0. 

Several  "housekeeping"  arrays  are  used  to  support  the  basis  data 
structures.  An  array  D  is  maintained  to  store  the  depth  of  each  node. 
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Jensen  and  Barnes  [Ref.  1]  seem  to  dismiss  the  preorder  traversal  scheme 
of  basis  representation  for  generalized  networks  due  to  the  apparent 
difficulty  of  assigning  depths  to  the  nodes  of  a  cycle.  This  difficulty 
is  overcome  by  defining  the  depth  of  all  cycle  nodes  to  be  0.  The  depth 
of  a  node  not  on  a  cycle  is  one  more  than  the  depth  of  its  predecessor. 

Additional  node-length  arrays  include  U  which  contains  the  values  of 
the  dual  variables  (or  simplex  multipliers)  and  X  which  contains  the 
value  of  the  right-hand  side  for  each  node's  conservation  of  flow  equation. 

A  node-length  array  FAC  is  used  to  store  cycle  factors  and  array 
PC  stores  information  about  the  incoming  non-basic  column  and  is  used  to 
effect  a  simplex  pivot.  These  arrays  will  be  discussed  later. 

To  associate  arc  numbers  with  the  node  pairs  connected  by  an  arc, 
we  maintain  two  arrays:  T  and  H.  T  is  an  arc-length  array  which  stores 
the  tail,  or  source  node,  of  each  arc.   If  the  arcs  in  T  are  sorted  so 
that  all  arcs  with  the  same  head  node  are  listed  contiguously,  a  space 
savings  can  be  effected  by  only  storing,  for  each  node  i,  the  location 
of  the  first  arc  whose  head,  or  destination,  is  i.  This  information  is 
stored  in  a  node-length  array  H;  H(i)  contains  the  index  of  the  first  arc 
in  T  whose  head  is  node  i.  Thus,  the  tails  of  all  the  arcs  whose  head  is 
node  i  are:  T(H(i)),  ...  T(H(i  +  1)  -  1).   If  H( i )  =  H( i  +  1) ,  then  no 
arcs  terminate  at  node  i. 

Associated  with  each  arc  are  arrays  which  describe  the  arc's  character- 
istics. These  arrays  are  naturally  indexed  by  arc  number.  The  array 
CP  stores  the  capacity  or  upper  bound  on  flow,  C  contains  the  cost  per 
unit  flow,  and  MUL  contains  the  arc  multiplier  or  gain. 
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It  is  only  necessary  to  consider  upper  bounds  in  the  solution  pro- 
cedure as  long  as  we  "transform  out"  any  lower  bounds  required  by  the 
model.  For  an  arc  k,  from  node  i  to  node  j  with  lower  bound  lb.  ,  this 
transformation  is  easily  performed  as  follows: 

CP(k)  <==  CP(k)  -  Ib,^ 

Obj  <==  Obj  +  Ibj^  *  C(k)      where  Obj  is  the  value  of 

the  objective  function; 

if  i  /  j 

X(i)  <==  x(i)  -  Ib^ 

X(j)  <==  X(j)  +  Ib^  *  MUL(k); 
if  i  =  j 

X(i)  <==  x(i)  -  Ib^  *  MUL(k). 

An  arc  k,  out  of  the  basis  at  its  upper  bound,  is  "reflected"  by  logically 
replacing  its  flow  variable  x.  by  (cpj^^  -  x,^)  and  that  reflection  is 
recorded  using  the  sign  bit  of  CP,  reminiscent  of  bounded  variable  simplex 
methods  [e.g.,  Ref.  3],  Figure  9  summarizes  the  suggested  arrays. 
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III.  ALGORITHM  DESCRIPTION 

As  with  any  simplex-based  algorithm,  the  solution  of  the  generalized 
network  problem,  (GNP),  involves  three  fundamental  operations:  priceout, 
ratio  test,  and  pivot.  Non-basic  variables  (arcs)  are  individually 
examined  in  the  priceout  to  determine  if  inclusion  in  the  basis  can  yield 
a  better  value  of  the  objective  function.  Given  a  favorable  incoming 
variable,  a  ratio  test  is  performed  to  determine  whether  changing  ,flow  on 
the  incoming  variable  will  result  in  the  incoming  variable  achieving  its 
opposite  bound  or  a  basic  variable  being  driven  to  one  of  its  bounds. 
The  variable  first  reaching  a  bound  is  deemed  "outgoing."  Finally,  the 
incoming  variable  replaces  the  outgoing  variable  in  the  basis  via  the 
pivot  operation.  We  examine  each  of  these  fundamental  steps  and  their 
specialization  to  generalized  networks  using  the  representations  and  data 
structures  presented  in  the  previous  section. 

First,  however,  it  is  appropriate  to  review  some  of  the  principles 
of  the  revised  simplex  method  [e.g.,  Refs.  16,  17].  Consider  the 
following  LP: 

(LP' )     min  ex 

s.t.  Ax=b       (Aismxn) 
0  <_  X  £  cp. 

The  matrix  of  technological  coefficients  A  may  be  partitioned  into 
a  basic  square  sub-matrix  B  (i.e.,  B  is  a  basis)  and  a  non-basic  sub- 
matrix  N.  If  B  consists  of  the  first  m   columns  of  A,  then  A  =  [B,  N]. 
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Similarly,  we  may  partition  the  variable  and  cost  vectors  x  and  c  as 
X  =  (Xg,  Xj^)  and  c  =  (Cg,  Cj^) .  (LP')  becomes: 


mm  ex  =  CgXg  +  c^x^ 


s.t.  Ax  =  [B,  N]rxg1=  Bxg  +  Nx,^  =  b 
x  >  0. 


The  upper  bounds  in  the  original  formulation  may  be  taken  into 
account  by  "reflecting"  any  variable  at  its  upper  bound.  That  is,  the 
value  of  a  reflected  variable  is  measured  as  the  "distance"  from  its 
upper  bound  cp.  (as  opposed  to  the  more  natural  way  of  measuring  values 
as  distance  from  the  lower  bound  0).  Any  variable  at  its  upper  bound 
cp,  is  replaced  by  (cpj^  -  Xj^);  to  effect  this  change,  coefficient 

column  k  is  replaced  by  its  negative  -A  and  the  right-hand  side  is 

i< 
transformed  to  b  -  A  cp.  . 

The  solution  to  the  LP  is  obtained  as 


BXg  +  Nxj^  =  b 

xg  =  3-^b  -  B-^Nx^. 


Substituting  this  into  the  cost  equation: 


=  CgB-^b  -  CgB-^Nx^  +  c^x^ 

=  CgB'  b  +  (Cj^  -C5B"  N)Xj^   (the  cost  of  the  current 

solution  in  terms  of  Xj^)  . 

Every  basic  solution  to  the  LP  has  the  characteristic  that  x..  =  0 
and  thus  the  cost  of  that  basic  solution  is  CoB~  b.  By  changing 
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1) 


the  current  basic  solution  and  thereby  making  an  element  of  x  non-zero, 

the  cost  of  the  solution  would  change  by  the  amount  (Cj^  -  CnB~  N)X|^. 

If  Cm  -  CoB'  N  <  0,  then  the  overall  cost  would  be  reduced.  Thus, 

elements  of  Xj^  for  which  C|^  -  CoB'  N  <  0  are  desirable  candidates 

to  enter  the  basis,  c.,  -  CqB"  N  are  called  the  reduced  costs.  CnB' 

N    B  B 

is  called  the  dual  solution  or  simplex  multiplier  set. 

Reviewing  the  algebraic  development  of  the  ratio  test,  we  begin 
with  the  general  expression  for  the  basic  variables  (1). 

Xg  =  B"-^b  -  B"^Nx^ 
If  X.  £X|M  is  the  incoming  non-basic  variable,  we  obtain 

Xg  =  B"-^b  -  z'^Xj^;  where  Z*^  =  B'^n'^  is  the 

current  incoming  column. 

To  preserve  feasibility,  each  of  the  x's  must  remain  within  its  respective 

bounds  (0,  cp,).  The  three  constraining  conditions  searched  for  by 

the  ratio  test  are  that  the  incoming  variable: 

(a)  reaches  its  opposite  bound: 

(b)  drives  a  basic  variable  to  its  opposite  bound: 
'b.  -  ^ik^k  =  ^PS.      ^^"^  ^ik  <  Q 

\   =  ^^3.  -  ^PB.)/^ik 

(c)  drives  a  basic  variable  to  its  current  lower  bound: 

^B.  -  ^ik\  =  0     ^°^  ^ik  >  0 

^k  ~  ^B^^ik       where  B^-  indicates  the  variable  which  is  basic 

in  row  i  of  the  constraint  matrix  and  z.,  is  the 

th  k 

corresponding  i   element  of  Z  . 
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The  pivot  operation  updates  the  current  solution  to  exhibit  the 
exchange  of  basis  elements.  This  is  accomplished  by  updating  the  represen- 
tation of  B"  and  the  right-hand  side  B"  b  via  a  pivot,  or  elementary 
transformation  operation  (as  B"  is  seldom  explicitly  extant),  using 
the  incoming  column  Z  .  If  the  incoming  variable  reaches  its  opposite 
bound,  then  no  basis  exchange  is  required;  however,  that  variable  must  be 
reflected. 

Tableau  arithmetic  is  not  the  most  efficient  method  for  manipulating 
the  problem  representation  of  a  linear  program,  nor  is  it  as  insightful 
as  the  matrix  arithmetic  of  the  revised  simplex  method.  As  shown,  the 
current  inverse  of  the  basis,  B"  ,  carries  with  it  enough  information 
to  generate  any  current  solution.  The  fundamental  principle  of  the 
revised  simplex  method  is  that  using  only  the  current  B"  and  the 
original  problem  coefficients  is  much  more  efficient  than  manipulating  a 
complete  tableau. 

A  further  enhancement  to  this  method  is  made  by  not  storing  B" 
explicitly  but  by  generating  it  dynamically  as  the  product  of  elementary 
vectors  (which  may  be  efficiently  stored  due  to  their  sparse  nature). 
These  same  principles  will  be  adhered  to  in  the  following  specialization 
of  the  (revised)  simplex  method  to  generalized  networks.  The  mechanism 
to  update  our  solution  will  be  based  on  the  dynamic  generation  of  the 
equivalent  to  a  revised  simplex  B~  . 

A.  PRICEOUT 

As  expected,  the  pricing  of  non-basic  variables  is  simply  the 
specialization  and  simplification  of  simplex  pricing.  The  reduced  costs 
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of  the  non-basic  variables  are  c  -  CqB"  N.  Let  u  =  CpB"  (the  dual 
solution).  Then  the  reduced  cost  (re)  for  a  non-basic  arc  k  is: 

N  is  a  column  of  the  network  constraint  matrix  having  at  most  two  non- 
zero entries,  i.e.,  +1  (in  the  row  corresponding  to  the  tail  of  arc  k) 
and  -m.  (corresponding  to  the  head  of  arc  k). 

The  reduced  cost  of  arc  k  oriented  from  i  to  j  simplifies  to 

"S  =  ^k  -  ^-  ^Vy 

or  in  the  data  array  notation  of  Section  II 

rc,^  =  C(k)  -  U(i)  +  MUL(K)  *  U(j). 

If  arc  k   is  a  reflected  arc   (denoted  by  CP(k)   <  0),  then  the  sign  of 

k    • 
re,    is  reversed.     (In  reflecting   a  variable  Xj^,  column  N     is 

multiplied  by  -1,   the  proper  cost   associated  with   that  variable   is, 

therefore,   -C(k)    and  the  non-zero  coefficients    in  column   N     change  to  -1 

and  +m.  ) . 

The  pricing  simplification  is  one  of  the  key  computational  advantages 
of  network  models.  At  most  one  multiplication,  one  addition,  and  one 
subtraction  is  required  to  price  a  non-basic  variable.   In  the  case  of 
pure  networks  m.  =  1,  the  multiplication  may  be  avoided. 

As  discussed  in  Section  II,  the  components  of  the  generalized  network 
basis  are  triangular  or  nearly  triangular.  As  a  consequence,  the  value 
of  the  dual  variables  may  be  found  by  forward  substitution  in  the  system 
of  equations  uB  =  Cn .  This  is  a  simple  matter  for  triangular  elements 
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of  the  basis  and  once  one  variable  is  determined  an  equally  simple 
matter  for  nearly  triangular  basis  elements.  A  detailed  discussion  of 
the  computation  of  dual  variables  will  be  presented  with  the  pivot 
operations. 

It  is  common  practice,  in  didactic  presentations  of  the  simplex 
method,  to  select  the  variable  with  the  most  negative  price  as  the 
incoming  basis  element.  This  approach  requires  pricing  all  the  non-basic 
variables  and  can  be  quite  time-consuming,  especially  in  network  models 
where  the  number  of  arcs  far  exceed  the  number  of  nodes.  In  fact,  any 
variable  which  prices  favorably  will  suffice  as  the  incoming  candidate 
and  any  pricing  mechanism  which  guarantees  discovery  of  all  favorably 
pricing  non-basic  variables  will  operate  correctly. 

Several  pricing  strategies  have  been  suggested.  It  has  been  shown 
that  selection  of  the  best  (most  negative  price)  from  a  limited  number  of 
favorably  priced  candidates  is  superior  to  simply  choosing  the  first 
favorably  priced  arc  [e.g.,  Ref.  18  for  LP  and  Ref.  5  for  GNP]. 

Some  popular  pricing  strategies  are:  (a)  price  out  the  first  g 
candidates  where  g   is  less  than  the  number  of  non-basic  variables;  (b) 
maintain  a  candidate  list  as  adopted  by  Glover,  et  al  .  [Ref.  5]  and  by 
Mulvey  [Ref.  19];  and  (c)  maintain  a  candidate  queue,  a  strategy  which  is 
employed  by  Bradley,  Brown,  and  Graves  [Ref.  6]. 

The  candidate  list  strategy  maintains  a  list  containing  at  most 
Q\   candidates.  Each  candidate  is  the  most  negative  arc  originating  from 
a  node.  Every  g2   pivots  (or  when  no  candidate  prices  favorably)  the  list 
is  refreshed.   If  less  than  gl  favorable  candidates  are  found,  then  gl  is 
set  to  the  number  of  favorable  candidates  and  g2   is  halved  (unless  gZ   =  1' 
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Mulvey  presents  computational  results  to  aid  in  the  selection  of 
gl   and  g2. 

The  candidate  queue  is  a  dynamic  list  of  "interesting  nodes"  and 
"good  arcs."  The  queue,  as  described  by  Bradley  et  al .  [Ref.  6],  employs 
two  arrays,  NSA  which  indicates  the  head  node  and  ISA  which  indicates  the 
tail  node  of  candidate  arcs.  ISA(k)  =  0  and  NSA(I<)  =  j  indicates  that 
node  j  is  an  "interesting"  node.  Arc  entries  are  derived  from  a  scan  of 
all  arcs  incident  to  a  node,  from  which  the  most  negatively  priced  arc  is 
entered.  The  candidate  queue  is  initialized  by  placing  demand  or  sink 
nodes  (right-hand  side  <  0)  on  the  queue  as  interesting  nodes.  (If  no 
such  nodes  are  found,  the  queue  is  initialized  with  any  node.)  A  general 
scan  of  a  head  node  will  select  the  most  favorably  priced  incident  arc 
and  place  that  arc  on  the  candidate  queue.  Initially,  favorable  prices 
are  most  commonly  caused  by  inherent  infeasibil ity  (i.e.,  paths  do  not 
yet  exist  between  supply  and  demand  nodes).  For  the  first  nns   pivots, 
known  as  the  opening  gambit,  the  head  and  tail  of  each  incoming  basic  arc 
may  be  entered  on  the  queue  as  "interesting"  nodes. 

For  each  pivot,  the  incoming  candidate  is  determined  by  oricinq  out 
nne   candidates  (including  all  of  the  arcs  incoming  to  an  interesting 
node).   If  no  favorable  arc  is  found  within  the  rme   candidates  examined, 
another  nne   are  examined,  and  so  forth.  If  the  queue  is  emptied  during 
the  first  nns   pivots,  the  queue  is  refreshed  by  a  general  scan  of  ipg   (a 
page)  of  head  nodes. 

After  the  end  of  the  opening  gambit,  the  queue  is  refreshed  after 
each  cycle  (pass)  through  the  queue  by  scanning  another  page  of  the  head 
nodes.  As  each  arc  is  priced,  the  favorably  priced  candidates  are 
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retained  in  the  queue,  interesting  nodes  are  replaced  by  their  most 
favorable  incoming  arc  and  unfavorable  candidates  are  dropped  from  the 
queue.  This  pricing  strategy  finally  concludes  by  pricing  every  non- 
basic  arc  (as  eyery   pricing  strategy  must)  to  ensure  terminal  conditions 
are  met  for  optimal ity.  Suggested  values  for  ipg,  nne,   and  nns   ^tq   given 
by  Bradley,  Brown,  and  Graves  [Ref.  6]. 

The  results  of  any  pricing  strategy  will  be  the  identification  of  the 
head,  tail,  and  arc  index  of  a  favorably  priced  incoming  arc  (or  the 
determination  that  there  are  no  favorably  priced  non-basic  arcs). 

B.  RATIO  TEST 

At  this  point,  an  incoming  arc  has  been  identified  and  the  algorithm 
must  now  select  an  appropriate  outgoing  member  of  the  basis.  Conceptually 
the  ratio  test  identifies  the  first  arc  whose  flow  would  reach  one  of  its 
bounds  as  flow  is  incremented  on  the  incoming  arc.  The  ratio  test  thus 
identifies  an  arc  whose  flow  would  either  increase  to  its  upper  bound  or 
decrease  to  zero  if  the  fullest  advantage  were  taken  of  the  favorable 
price  of  the  incoming  arc.  The  incoming  arc  may  well  be  the  constraining 
arc  identified  by  the  ratio  test  and  therefore,  the  basis  would  not 
change;  however,  network  flows  would  be  updated  with  the  flow  on  the 
incoming  arc  being  "sent"  to  its  opposite  bound  by  reflection. 

For  an  incoming  arc  k,  the  ratio  test  searches  for  the  minimum  of: 

(a)  cp,^ 

(b)  (xg_  -  cPb_)/z^,^      for  z^^  <  0    V  Xg^Xg 

(c)  X3_/z.^  for  z.^  >  0    V  Xg^eXg^ 
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own . 


We  know  cp.  j  =  1,  ...  ,  n  and  the  current  basic  solution  Xn .  Therefore, 

the  ratio  test  can  easily  be  performed  if  z^-.  i  =  1,  ...  ,  m   is  kn 

To  determine  Z  (the  current  incoming  column),  we  must  solve  the  system 

of  equations  BZ  =  N  (remembering  to  use  -N  if  arc  k  is  reflected). 

k    -Ik 
The  obvious  algebraic  solution  of  this  system  is  Z  =3  N  . 

We  can  generate  N  on  demand  since  it  is  the  coefficient  column  associated 

with  arc  k  (oriented  from  node  i  to  node  j).  We  determine  i  and  j  from  the 

H  and  T  arrays  and  thereby  determine  the  two  non-zero  entries  of  N  as  +1 

(corresponding  to  i)  and  -m.  or  MUL(k)  (corresponding  to  j). 

We  do  not  have  an  explicit  representation  of  B'  and  thus  must 

resort  to  indirect  methods  of  solution.  It  is  in  this  endeavor  that  the 

predecessor  relationship  proves  to  be  extremely  convenient.  Using  the 

predecessor  array,  it  is  possible  to  solve  for  Z  directly  by  substitution 

k    k 
in  the  triangular/nearly  triangular  system  of  equations  BZ  =  N  . 

N  is  a  non-basic  column  of  the  constraint  matrix  and  therefore 

has  at  most  two  non-zero  entires:  +1  and  -m.  .  N  can  be  expressed 

as: 

n''  =  e.  +  (-m^ej) 

where  e-  and  e-  are  the  i   and  i   unit  vectors  corresponding 
to  the  non-zero  entries  of  N  .  Solving: 


8QJ  =  -Vj. 


will  lead  us  to 
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=  BQ^  +  BQ'^' 

=  B(Q^  +  qJ) 
thus,  7*^  =  0^+  Q"^. 
The  entering  arc  connects  nodes  i  and  j.  Define  the  "i-backpath"  as 
the  path  between  node  i  and  its  root/cycle.  Similarly,  define  the 
"j-backpath"  as  the  path  between  node  j  and  its  root/cycle.  It  will  be 
shown  that  Q^  and  Q^^  correspond  conceptually  to  the  i-backpath  and 
j-backpath,  respectively. 

Remember,  the  basis  is  made  up  of  disjoint  components  each  of  which 
is  triangular  or  nearly  triangular.  Working  with  triangular  components 
first,  we  find  for: 

BQ  =  e.    (or  similarly,  -ni.e.) 
that  elements  of  Q  beyond  the  i   component  must  be  zero.  This  can  be 
seen  by  solving  the  triangular  system.  If  e-  is  of  dimension  h  x  1  and 
i  <  h,  then  e.(h)  =  0,  thus  Q^h)  =0.  If  i  <  h  -  1,  then  e.(h  -  1)  =  0, 
and  since  Q^(h)  =  0,  then  Q^(h  -  1)  =  0.  The  same  reasoning  applies  to 
Q^i  +  1)  =  Q^i  +  2)  =  ...  =  g^h)  =  0.   The  i^^  element  of  0^  must  be  +1 
or  -m.  .   If  there  is  an  entry  in  the  triangular  matrix  for  the  oredecessor 
of  i,  which  is  in  the  i   column,  it  will  oroduce  a  non-zero  multiplica- 
tion and  thus  must  be  offset  by  an  entry  in  the  P(i)  row  of   q\  since  e. 
has  only  one  non-zero  entry.  That  "offsetting"  element  may  also  have  a 
predecessor  and  thus  another  non-zero  entry  in  Q  will  be  appropriate. 
Q^  will  therefore  have  non-zero  entries  in  rows  i,  P(i),  P(P(i)),  etc., 
(i.e.,  the  i-backpath).  The  preceding  argument  also  applies  to  Q"^ , 
which  will  have  non-zero  entries  in  only  rows  j,  P(j),  etc. 
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To  see  this  and  to  determine  the  non-zero  entries  of  Q,  consider 
the  triangular  system  in  Figure  10. 

^root 


n 


I 

0 

0 
0 


P(P(i)) 


P(il 


'P(P(P(i))) 


ap(p(i))    ^?(P(i)) 


^P(i)    ^P(i) 


Figure  10.  Triangular  Basis  Component 


The  values  of  a  and  b  are  dependent  on  the  original  arc  orientation 
(recorded  by  the  sign  bit  of  the  predecessor  array  P;  if  P(j)  >  0,  then 
a  =  +1  and  b  is  the  negative  of  the  fnultiolier  value  (-m.  );  if 
P(j)  <  0,  then  a  =  -m,  and  b  =  ^i,  '-/  j  =  i ,  ^(  i) ,  ...  ,  root.   Solving: 


i-1 


a- 


i-1 


i   ! 


'i-1 


0   1 


mere   s  =  1  or  -m.  , 
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we  obtain: 


and  in  general 


q,-  =  s/a^. 

%   =  -  ^VlVl^/^h'    0  <  h  <  i. 


If  we  are  dealing  with  a  nearly  triangular  basis  element,  we  again 
find  that  the  only  non-zero  entries  in  Q^  occur  in  entries  corresponding 
to  the  backpath  between  node  i  and  the  "root"  cycle.  Suppose  the  entering 
arc  is  incident  to  node  i  which  is  on  the  cycle  shown  in  Figure  11. 


Figure  11.  Cycle 

The  backpath  is  thus  i,  P(i),  P(P(i))  ...  r  where  P{r)   =  i.  The  nearly 
triangular  basis  component  is  as  follows. 


'P(P(?( 
^P(P(i) 


'P(P(i)) 


7(i 


'5^(1) 


'P(r) 
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Solving: 


^0 

^1 


b. 


i-1 


a,- 


i-1 


•^i-l 


(2) 


we  obtain 


an 


where  f  = 


where  s  =  1  or  -m,  and  the  a's  and 
b's  are   as  before; 


q^  =  (s/a.)(l/f) 

^  %  =  -  (Vi%+i)/^h, 


0  <  h  <  i 


1  -   n   (b^/a^^). 
h  =  0 


The  solution  procedure  is  described  by  Oantziq  [Ref.  3].  The  nearly 
triangular  system  can  be  easily  solved  with  one  of  the  variables  aopearinq 
as  a  solution  parameter  thereby  creating  a  triangular  system.  That 
parameter  can  then  be  determined  as  the  solution  to  a  pair  of  two-variable 
simultaneous  equations.  Kennington  and  Helgason  [Ref.  2]  give  a  detailed 
derivation  of  the  solution.  A  similar  solution  technique  is  illustrated 
herein  with  the  discussion  of  the  dual  update. 
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These  equations  are  the  same  as  in  the  triangular  case  except  for  the 
use,  once,  of  the  cycle  factor  f .  The  cycle  factor  (or  loop  factor, 
Dantzig)  is  the  same  for  every  node  on  the  cycle  and  may  be  computed  when 
the  cycle  is  created.  The  array  FAC  is  used  to  store  the  cycle  factor  so 
that  it  is  available  for  immediate  use.  The  update  of  FAC  will  be 
discussed  in  the  pivot  section.  Again,  the  above  solution  procedure  for 
Q^  applies  equally  well  to  determining  Q"^ . 

We  have  approached  this  problem  as  if  there  are  two  distinct  backpaths. 

However,  the  incoming  arc  may  be  such  that  the  two  backpaths  converge.  The 

node  at  which  the  backpaths  meet  is  known  as  the  "join."  The  two  separate 

systems  of  equations  must  be  solved  up  to  the  join  using  the  iterative 

method  discussed  above.  The  separate  solutions  for  q.  .  obtained  from 

^  join 

the  converging  backpaths  should  be  added:  q.  .  =  q-  ■   .  +  q.  . 

^  ^     ^  join   ^join,  1   ^join,  j 

and  only  one  consolidated  backpath  need  be  computed  after  the  join 
using  the  same  interative  formula: 

-^h^l%^l 

%  = r-  ' 

h 

The  q's  obtained  are  the  elements  of  Z  we  need  to  compute  ratios. 

For  use  in  the  pivot,  we  store  these  values  in  the  array  PC.  The  array 

-1  k 
PC  contains  B  N  ,  which,  in  revised  simplex  terms,  is  the  incoming 

non-basic  column  updated  by  the  current  B~  . 

Obviously,  detection  of  a  join  is  extremely  important  in  reducing 

the  number  of  computations  required  to  compute  ratios.  The  depth  array  D 

is  used  to  help  detect  a  join  in  the  following  way.  Let  arc  k,  joining 

node  i  and  node  j,  be  the  entering  arc;  let  i  refer  to  the  node  which  is 
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currently  being  "visited"  on  the  i-backpath  and  let  j  indicate  the  node 
currently  "visited"  on  the  j-backpath. 

We  begin  with  i  and  j  as  the  origin  and  terminus  of  arc  k.  The 
depths  of  i  and  j  are  compared  and  the  backpath  of  the  one  with  the 
greater  depth  is  iterated  using  the  predecessor  relationship.  Ratios  and 
the  non-zero  elements  of  Z  are  found  one  at  a  time  as  each  node  on  the 
backpath  is  visited.  If  the  j-backpath  is  deeper,  we  iterate  it  back 
until  the  depth  is  the  same  as  the  i-node.  The  current  i-  and  j-nodes 
are  checked  to  see  if  they  are  the  same  node.  If  they  are,  the  current 
element  of  Z  is  computed  as  the  join  and  the  common  backpath  is 
iterated,  continuing  to  search  for  the  minimum  ratio.   If  the  two  nodes 
are  not  the  same,  then  the  i-backpath  and  the  j-backpath  are  both  iterated 
once  and  a  join  is  checked  for  again. 

In  summary,  the  mechanics  of  the  ratio  test  are:  if  arc  k(i,  j) 
is  the  entering  arc,  determine  the  minimum  among: 

1.  cp^ 

2.  (Xg_  -  cpg_)/z.^       for  z.^  <  0    V  Xg^cXg 


1 


3.  Xg/z.^  for  z.^  >  0    V  Xg.Xg. 

1  1 

To  do  this: 

(a)  Determine  cp,  as  C?(k). 

(b)  Determine  whether  the  i-backpath  or  j-backpath  is  deeper  by 
checking  the  depths  of  nodes  i  and  j,  the  origin  and  terminus  of  the 
entering  arc.  Let  i  refer  to  the  node  currently  being  visited  on  the 
i-backpath  and  j  refer   to  the  current  node  on  the  j-backpath.  Using  the 
predecessor  relationship,  iterate  up  each  backpath,  oerforming  the  steps 
described  below  on  each  node  encountered.  When  the  current  i  and  j 
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nodes  are  of  the  same  depth,  check  for  a  join.   If  a  join  is  found,  add 

the  values  of  "q"  for  that  node  obtained  from  each  separate  backpath  to 

obtain  q.  •  .   If  a  join  is  not  found,  iterate  up  one  level  on  both  the 
join       ^  '^ 

i-  and  j-backpath  and  check  for  a  join  again. 

(c)  Determine  the  right-hand  side  of  the  triangular/nearly  triangular 
system  of  equations  (2).  This  is  accomplished  by  setting  s  =  1  for  the 
i-backpath  (starting  at  the  tail  of  the  entering  arc)  and  setting  s  =  -m. 
for  the  j-backpath  (starting  at  the  head  of  the  entering  arc)  in 
equation  (2) . 

(d)  For  each  node  on  the  i-  and  j-backpath,  determine  the  sign  of 
the  current  node's  predecessor  P(i).  Use  this  to  determine  a  and  b  in 
equation  (2)  as  P(i)  >  0  =>  a  =  1,  b  =  -my^;  P(i)  <  0  =>  a  =  -m|^,  b  =  1 
where  m,  is  the  multiplier  of  the  arc  joining  i  and  P(i). 

(e)  Compute 

q-  =  s/a^-  if  i  is  the  first  node  on  the  backpath; 
otherwise, 

q.  =  -  (b.^^q.^^)/a.,  where  P(i+1)  =  i. 

(f)  Determine  if  a  join  exists;  if  so,  add  the  i-  and  j-backpath 

values  of  q  to  obtain  q-  •  . 

join 

(g)  Determine  if  the  current  node  is  the  first  node  encountered 
on  the  root  cycle  (D(i)  =  0).   If  so,  multiply  q-  by  l/FAC(i}. 

(h)  The  q-  thus  determined  is  the  z-j^  required  to  perform  the 
ratio  test 

^B/^ik  '^   ^ik  >  0 

(xg_  -  cp3_)/Zi^       if  z^^  <  0 

cpg  is  found  in  CP(IVAR(i)) . 
i 

43 


(i)  The  ratio  test  may  be  terminated  when  a  zero  ratio  is  found  (a 
de  facto  winner)  or  we  completely  iterate  the  i-  and  j-backpath  performing 
ratio  tests.  The  end  of  a  backpath  is  signified  when  the  first  node  with 
depth  zero  is  encountered  for  the  second  time.  (For  a  rooted  tree,  the 
root  is  its  own  predecessor  and  will  be  "encountered"  twice  in  immediate 
succession  by  following  the  predecessor  relationship.  On  the  other  hand, 
the  nodes  of  a  cycle  all  have  depth  zero  and  once  the  cycle  is  completely 
iterated,  you  will  find  a  previously  encountered  node  of  depth  zero.) 

C.  PIVOT 

We  have  now  determined  the  entering  arc  (priceout),  the  leaving 
arc  (ratio  test),  and  have  generated  the  entering  column  and  stored  it  in 
array  PC.  The  basis  representation  and  arc  flows  must  now  be  updated. 

If  CP(k)  is  the  minimum  ratio,  then  the  incoming  arc  has  been  deter- 
mined to  be  more  constraining  than  any  element  of  the  basis.  Therefore, 
the  incoming  arc  k  is  also  the  outgoing  arc.  The  update  is  accomplished 
by  reflecting  arc  k.  The  original  basis  remains  unchanged  and  only  the 
right-hand  side  need  be  updated.  The  general  expression  for  the  right-hand 
side  was  derived  in  equation  (1) 

Xg  =  B'^b  -  B"^Nx^, 

where:  3"  b   is   the  old  right-hand  side 

Xj^       are  the  values  of  the  non-basic  variables. 

The  non-basic  variables   are   all    zero  except,   conceptually,   the   incom- 
ing variable  which    is   set   at    its   upper  bound  CP(k).     Therefore,   only  the 
k       column   of  B"   N"^    is   required,   which    is   precisely  the   information 
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now  stored   in  PC  as   an   artifact  of  the  ratio  test.     The  right-hand 
side  (rhs)   update   is  then: 


rhs.   =  old  rhs-   -  PC(i)   *  CP(I<) 

'  '  (Vi) 

X(i)   <==  x(i)   -  PC(i)   *  CP(k). 


If  the  ratio  test  determines  that  an  exchange  of  basis  elements 
is  required,  a  more  involved  update  procedure  takes  place.  The  basis 
representation  found  in  P  and  IT  must  be  updated,  as  well  as  the  dual 
variables  U  and  the  flow  X.  The  FAC  array  must  be  maintained  for  nodes 
on  cycles  and  the  IVAR  array  contains  the  index  of  the  basic  arc  associ- 
ated with  each  node.  As  such,  IVAR  must  also  be  updated  during  the 
pivot. 

To  conceptually  understand  the  basis  update  procedure,  return  to 
the  graphical  representation  of  the  basis.  The  root/cycle  of  each 
basis  component  is  drawn  at  the  top  and  the  immediate  successors  of  each 
node  are  depicted  below  that  node  (as  with  a  genealogical,  vice  a 
biological ,  tree) . 

It  is  important  to  carefully  distinguish  between  "types"  of  successors 
The  "successors"  of  a  node  will  mean  those  nodes  immediately  encountered 
on  all  paths  from  a  node  except  the  backpath.  "Preorder-successors"  will 
indicate  nodes  determined  t^y   iteration  of  the  preorder  function  IT,  Arcs 
are  drawn  to  indicate  predecessor  relationships  (i.e.,  the  arcs  always 
point  "up"),  with  the  sign  of  the  predecessor  being  used  to  record  the 
"correct"  orientation  of  the  arc. 
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Figure  12.     Pre-Pivot  Generalized  Network  Basis 

The  simple  basis  pictured  above  will   be  used  to   illustrate  the  basis 
update  mechanism  when  the   incoming  arc   is  not  the  "winner"  of  the  ratio 
test.     Let  the  incoming  arc  be  design-ated  as  '<r(i,  j)   and  the  outgoing 
arc  be  k,(v,  w) .     Assume  the  outgoing   arc   is  on  the  j-backpath   (if   it 
is  not,  reorient   arc  k  so  that   it   is).     Similarly,   let  node  v  precede 
node  w  on  the  j-backpath.     Assume  the  entering  arc   is  kc(2,   8)    and  the 
leaving  arc   is  k, (5,   6).     The  backpath  from  the  destination  node  j   of 
the   incoming  arc  to  node  v  of  the  outgoing  arc   is  called  the  "j-stem." 
The  j-stem  in  the  example   is   along  nodes  3,   6,   7,   5.     Dropping  arc  k   (5,   6) 
and   adding   arc  kr-(2,   3)    forms   a  new  basis   as   shown    in  Figure  13. 


®^; 


aj) 


Figure  13.  3asis  Update  Example  1 
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Another  example  of  a  basis  update  is  shown  in  Figure  14. 


k      (6,  3) 

^^_  (2,  I) 

j  -  stem  3,  2 


Figure  14.     Basis  Update  Example  2 

The  update  of  any  basis  can  be  similarly  accomplished  by  application 
of  the  following  general   rules: 

1.  Reverse  the  arc  orientation  of  all    arcs  on  the  j-stem  (which 
lies  on  the  backpath  containing  the  leaving  arc). 

2.  Orient  the  entering  arc  such  that   it  precedes  the  j-stem  and 
has  the  same  orientation   as  the  redirected  j-stem. 

The  basis  update  shown   in  Figure  15  displays  how  a  new  cycle   is 

created. 


Figure  15.  Cycle  Creation 
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As  noted  in  the  ratio  test  section,  if  the  i-  and  j-backpaths  merge, 
the  node  at  which  they  merge  is  called  the  join.  If  the  leaving  arc  lies 
beyond  the  join,  then  a  new  cycle  is  formed,  as  is  the  case  in  the  above 
example.  The  new  cycle  becomes  the  root  of  its  basis  component  and 
that  component  tree  is  "rehung"  from  the  new  cycle.  The  backpath  from 
the  source  node  (i)  of  the  entering  arc  to  the  node  whose  predecessor  is 
the  join  is  known  as  the  i-stem.  In  this  example,  P(3)  =  5  =  join;  the 
i-stem  is  therefore  composed  of  nodes  4,  2,  and  3.  If  node  i  lies  on  the 
j-backpath,  the  i-stem  is  null. 

The  graphical  display  of  a  basis  update  must  be  translated  into  an 
algebraic  update  of  the  data  structures  discussed  in  Section  II.  Each 
node  affected  by  the  basis  update  must  be  "visited"  by  the  algorithm  and 
the  associated  array  values  modified.  Clearly,  it  is  advantageous  to 
visit  a  node  only  once.  A  one-pass  update  can  be  achieved  as  described 
below. 

The  pivot  algorithm  iterates  up  the  j-stem  node  by  node,  updating 
the  arrays:  X,  IT,  P,  U,  IVAR,  and  0.   If  a  join  is  encountered  (the 
existence  of  a  join  is  predetermined  by  the  ratio  test),  the  algorithm 
switches  to  the  i-stem  and  iterates  up  the  i-stem.  Upon  completion  of 
the  i-stem,  the  algorithm  returns  and  completes  iterating  the  j-stem. 
The  stems  are  iterated  by  using  the  predecessor  function  P. 

The  depth  0  and  the  dual  U  must  be  updated  for  all  of  the  stem  nodes 
visited  by  the  algorithm,  as  well  as  for  their  preorder-successors . 

It  is  also  convenient  to  modify  IT  as  these  nodes  dre   visited.  The 
preorder-successors  of  a  node  can  be  divided  into  two  classes:  the  left 
preorder-successors  and  the  right  preorder-successors.  The  left  preorder- 
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successors  are  found  by  iterating  the  preorder  list  IT  from  the  current 
stem  node  to  the  stem  node  whose  predecessor  is  the  current  stem  node. 
The  right  preorder-successors  of  a  node  are  any  unvisited  nodes  found  by 
further  iteration  of  IT  up  to  the  end  of  the  tree  whose  root  is  the 
current  stem  node.  The  end  of  this  current  tree  is  found  by  checking  the 
depth  of  each  preorder-successor  of  the  stem  node.  If  that  depth  is  less 
than  or  equal  to  that  of  the  current  stem  node,  then  the  start  of  a  new 
tree  is  identified. 

The  update  of  IT  is  relatively  easy  because  IT  changes  only  for  nodes 
on  the  stem  and  their  last  left  and  right  preorder-successors. 

0.  IT  (last  right  preorder-successor  of  i)  =  j 

1.  j-stem  update  of  IT 

a.  IT  (last  left  preorder-successor)   <==  first  right  preorder- 

successor 

b.  IT  (last  right  preorder-successor)  <==  P  (stem  node) 
If  there  is  no  last  right  preorder-successor,  then 

IT  (stem  node)  <==  P  (stem  node) 

2.  i-stem  update  of  IT 

a.  IT  (last  left  preorder-successor)   <==  first  right  preorder- 

successor 

b.  IT  (last  right  preorder-successor)  <==  node  whose  predecessor 

is  current  stem  node. 

Step  lb  differs  from  2b  because  the  predecessor  relationship  for  the 
j-stem  must  be  reversed. 

Some  of  the  nodes  encountered  in  the  preorder  may  belong  to  the 
i-stem  and  its  preorder-successors;  these  nodes  must  be  avoided  during 
the  j-stem  update,  using  the  "preorder  link."  We  process  the  part  of  the 
j-stem  below  the  join  and  then  process  the  i-stem  (if  there  is  a  join). 
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The  preorder-successor  of  the  last  i-stem  node  is  recorded  and  this  node 
is  termed  the  "preorder  link."  Having  finished  processing  the  i-stem,  we 
return  to  the  j-stem  (or  common  backpath  above  the  join)  and  continue 
processing  stem  nodes  and  their  preorder-successors.   If  an  i-stem  node 
is  encountered,  the  remainder  of  the  i-stem  and  their  preorder-successors 
are  avoided  by  immediately  skipping  to  the  preorder  link.  This  subtle 
processing  twist  significantly  reduces  the  time  spent  searching  for  the 
preorder-successors  of  the  stem  nodes. 

The  i-stem  must  be  visited  only  if  there  is  a  join.   If  there  is  no 
join,  then  the  j-stem  is  appended  to  the  i-backpath  by  the  entering  arc 
and  only  the  new  preorder-successors  of  i  (i.e.,  the  j-stem  and  their 
preorder  successors)  need  be  updated. 

The  method  we  have  established  for  the  update  of  the  preorder-successor 
relationship  ensures  that  all  the  preorder-successors  of  a  (cycle)  node 
are  encountered  before  the  next  (cycle)  node.  While  traversing  a  stem,  a 
stem  node  may  be  encountered  with  a  right  preorder-successor  that  has 
depth  zero  and  is  therefore  on  the  root  cycle  (e.g.,  Figure  16).  We 
:herefore  observe  that  given  that  the  leaving  arc  is  on  the  bacioath  of 
the  j-stam,  a  cycle  which  is  being  broken  will  always  be  encountered  as 
the  right  preorder-successor  of  a  stem  node. 


Figure  16.  Preorder-Successor  with  Depth  Zero 
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1/ 
From  the  development  of  the  computation  of  Z   ,    it  can  be  seen  that 

the  arc  flows  will   only  change  on  the  backpaths.     The  array  X  represents 

the  basic   arc  flows   and  must  be  changed  for  each   arc  on  the  j-stem.      If 

the  entering  arc  forms   a  cycle,   flow  changes  occur  on  the   i-stem  as  well 

as  the  j-stem.     The  update   is  obtained  from  equation  (1): 

Xr   =  B"-^b  -  B'^Nx^ 

^  ^  (Vi) 
^^  X(i)    <==  X(i)    -  PC(i)   *  RATIO. 


Only  one  element  of  Xj^  (i.e.,  x.  )   affects  the  update  of  the  right-hand 
side.     Clearly,    if  no  cycle   is  created  and  the  minimum  ratio   is  zero,   the 
right-hand  side   is  not  changed. 

The  dual   variables,   stored   in  the  array  U,   are  determined  as 

u  =  C3B-I 
C3  =  uB. 
Therefore,  Cj^  =  u^   -  mj^u-         V  k(i,  j)    e  basic  arcs 
u.   =  c^  +  m^u. 
Uj  =   (c^  -  u^)/(-m^). 

Enforcing  the  above  relationship  will   determine  the  dual    variables   as 
the  stems   are  traversed   in  turn.      If  a  cycle   is  not  created,  only  the 
duals  for  the  j-stem  nodes    and  their  preorder-successors   need  be  computed 
For  a  node   i    and  associated  basic   arc  k,   the  update  depends  on  the 
orientation  of  arc  k. 

U(i)    =  C(k)   +  MUL{k)    *  U(P(i))  for  k(i,   P(i));   P{i)    >  0 

and 

U(i)    =   (C(k)   -  U(-  P(i)))/{-  MUL(k))  for  k(P(i),    i);    P(i)    <  0. 
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The  preorder  traversal  scheme  (coupled  with  the  reversal  of  the 
j-stem  predecessor  relationships)  ensures  that  the  dual  of  a  predecessor 
node  is  known  prior  to  its  use  in  updating  the  dual  of  an  immediate 
preorder-successor . 

When  a  new  cycle  is  created,  the  i-  and  j-backpaths  terminate  at  a 
new  root  (the  newly  created  cycle).  The  duals  for  this  entire  new  basis 
component  must  be  recomputed.  To  do  this,  the  dual  variable  must  be 
determined  for  one  of  the  cycle  nodes.  Algebraically,  the  situation  is 
analogous  to  the  determination  of  Z  during  the  ratio  test.  Once  the 
dual  of  one  cycle  node  is  established,  the  dual  variables  for  the  remaining 
cycle  nodes  and  their  preorder-successors  may  be  determined. 

The  ratio  test  gives  ample  warning  that  a  new  cycle  will  be  formed. 
If  the  leaving  arc  k. (v,  w)  lies  above  the  join,  a  new  cycle  will  be 
created  as: 


«  < 

1 
) 

Figure  17.   Leaving  Arc  Above  the  Join 

The  new  cycle  is  formed  (disolaying  the  new  Dredecessor  function)  as 
snown  in  Figure  13.  The  nearly  triangular  system  used  to  compute  the 
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p(p(l))      p(l)       l  =  p(j) 


Figure  18.     Predecessor  Update  for  New  Cycle 
dual   variables  corresponding  to  this  newly  formed  basis  component    is 
uB  =  c  = 


(Uq.      . 


'   V 


o       1 

g-1 

g-1         q 


=    (Cn, 


•V 


Dantzig  [Ref.  3,  3.  423]  describes  'now  to  solve  systems  such  as 

th  i  s , 

"...  by  treating  one  variable  of  the  loop  as  the  parameter  and 
evaluating  all  others  in  terms  of  it  as  we  proceed  around  the  looo. 
Upon  completion  of  this  circuit  a  second  expression  for  the  parameter 
will  result,  and  by  equating  the  two  expressions  we  :nay  evaluate  it 
numerical ly." 

Assume  the  following  cycle  has  been  formed: 
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Figure  19.  Four  Node  Cycle 


producing  this  nearly  triangular  system: 


(U^,  U3,  U2,  u^) 


'1 


"  V  c^,  c^ »  C2 »  C  j^ ) 


Performing  the  matrix  multiplication   and  solving  for  u: 

U4  =  (C4  -  b4U^)/a^ 

U3  =  (C3  -  b3U^)/a3 

Cn       I       Up  =  (Cp  -  bpU3)/a2 

u-j^  =   (Ct    -  b^U2)/a^ 

Choosing  u^   as  the  parameter,   we  obtain  through  forward  substitution 


^3^4  *  ^3^3  =  ^3 


^2^3  *  ^2^2 


b,U2  +  aiU,    =  Ci 


Un   =  -=~  ;c,   -  bvJ2) 


M 


So, 


=  17^'!   -17  ^'2  -  ^2^  17  '"3  -  '^3^4)))) 

b  P 


bi  b^  b^ 

c.    -  T^  (C-,  -  ~  (c-,   -  T^  (c,   -  b,u^))) 


^3    "^^ 


a,   -4 


'4-|, 
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b         b 


^3"^       ^S^'a'i 


)) 


=i-rc  A 

^1    L^       '2 


^3^4 


^3^4 


u,)] 


a^  L  1         32  a 


b,bpC-       b,b_b_c-       b,b„b^b 


^2^3  ^2^3^4 


iWi    1 

^2^3^4       ^  J 


y        —    ——    —    — ^— — .       T   — — —    —    — ^— — — —   T    -^-^_^-^_    (J. 

1       a-i       a-i  din        a-]  aoa-j       a-i  aoa-ia^       a-i  apa^a^     J- 


^1^2^3^4/     ^ 


^_  ^2     ^^l^.^Wi 


a-i       a-i  a^ 


a-i  a^a-i       a-i  aoa-ja/i 


e: — 5: 


1 


V^ 


^1       ^2       ^3       ^4 


The  denominator   is  precisely  the  cycle  factor  f  previously  defined   as 

1  -    n    -  — 
d=l      ^d 


The  sign  alternation  implied  by  the  terms  can  be  observed  in  the 

a 


derivation  of  u 


1 


Thus 


^1  -^1 


C^  -  D^ 


V 


.3  -  D3  \  a^ 


^  f 


Letting  ji  =  —  , 
'+   a^ 


C3  -  b3u' 
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^2 


h  -  Vi 


„    _  'l  ■  ^"2 


we  obtain 


u,   = 


^1  -"i 


r 

•^2 

^3-^3 
^            ^3 

"i 

)| 

V 

'2 

J 

X   ^ 


C3    -    bjU' 

'^1  '  "1  V       "2     /  ,  1 

— ^-^— — ^-^— — — ^— ^— ^^^—     A  -^ 

^1  ^ 


Vi¥il,  1 

—————     A      -^ 

ai  f 


U-,    = 


ui  x^  . 


Looking  at  a  similar  triangular  system 


(u;|,  u^,  u^,  u^) 


^4     ^3 


^3     H 


h  \ 


~     ^  ^4 '  ^3'  ''9»  ^\' 


we  find  that  it  is  easily  solved  as 


11"  =  _1 

^   ^4 


u 


c  -  b  u" 
3    3  4 


3 


c  -  b  u" 


c  -  b  u" 
...    1    12 


The  solutions  for  the  u"  are  the  same  expressions  assumed  for  u'. 
The  values  for  u'  may  therefore  be  obtained  as  the  solution  of  a  similar 
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triangular  system;  the  dual  variables  associated  with  a  cycle  (and  hence, 
a  nearly  triangular  system)  may  be  obtained  from  u' . 

To  obtain  the  key  dual  value  u,,  we  must  therefore  determine  ui, 
uJ,,  Up,  u|,  and  finally  Ui.  In  general,  if  we  have  a  cycle  composed 
of  nodes  1,  ...  ,  g,  the  determination  of  the  dual  of  one  of  these 
nodes,  u, ,  is  found  by  modified  forward  substitution  as  illustrated 
above.  With: 

u^  =  u|  X  (1/f) 

g 

where     f  =  1  -  n  -  ^^r^^r"^ 

r=l 

and     u'  =  cg/ag 

u;  =  (c^  -  b^u;^^)/a^    Vr  =  g-1,  ...  ,  1; 

c  is  the  cost  of  the  arc  associated  with  node  g. 

Theoretically,  we  may  start  at  any  node  on  the  cycle,  iterate  around 

the  cycle  computing  u'  and  accumulating  the  terms  of  n  ;  in  practice 

a 

we  choose  to  start  at  the  terminus  of  the  entering  arc  (i.e.,  the  j-node), 
After  visiting  all  the  cycle  nodes,  the  cycle  factor  f  is  computed  and 
stored  in  FAC  for  each  node  on  the  cycle  and  the  key  dual  variable  is 
determined. 

The  dual  of  i  is  the  keystone.  From  the  dual  of  i,  we  can  compute 
the  dual  of  j  and  those  of  the  entire  j-stem,  as  well  as  those  of  the 
i-stem. 

As  stated,  we  want  to  start  at  the  j-node  and  proceed  around  endinq 
at  the  i-node.  This  traversal  cannot  be  conveniently  supported  by   either 
the  pre-pivot  or  post-pivot  predecessor  relationships.   It  is  necessary 
to  follow  the  pre-pivot  j-stem  predecessor  relationship  and  the  reverse 
of  the  i-stem  predecessor  relationship. 
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To  perform  this  intricate  maneuver,  the  cycle  factor  and  key  dual 
value  of  a  newly  created  cycle  is  computed  immediately  after  the  ratio 
test  and  before  the  basis  update.  First,  the  pre-pivot  j-stem  predecessor 
relationship  is  iterated  from  the  j-node  to  the  join,  computing  a  partial 
product  of  the  cycle  factor.  The  i-stem  is  then  iterated,  performing  no 
computations  other  than  storing  the  reverse  predecessor  relationship.  (A 
convenient  place  to  do  this  is  provided  in  the  array  U,  which  now  contains 
obsolete  dual  values  and  must  be  recomputed  subsequently.)  Once  at  the 
end  of  the  i-stem,  the  reverse  i-stem  path  may  be  followed  to  complete 
the  computation  of  the  cycle  factor  and  the  dual  value  of  the  i-node. 
After  this  newly  formed  cycle  has  been  traversed,  the  pivotal  update  can 
be  accomplished,  following  both  the  i-  and  j-stem  and  iterating  IT  to 
update  preorder-successors . 

The  update  of  the  depth  array  is  straightforward.  If  a  new  cycle  is 
formed,  the  nodes  on  the  new  cycle  are  assigned  depths  of  zero.  The 
depths  of  the  left  and  right  preorder-successors  of  the  i-  and  j-nodes 
are  updated  as:  0(s)  =  D(P(s))  +  1  (after  the  update  of  the  j-stem  pred- 
ecessor relationship).   If  no  new  cycle  is  formed,  then  0(j)  =  0(i)  +  1 
and  the  preorder-successors  of  the  j-stem  nodes  are  updated  accordingly. 
Again,  the  preorder  ensures  that  D(P(s))  is  available  when  it  is  time  to 
compute  0(s) . 

IVAR  contains  the  index  of  the  basic  arc  associated  with  each  node. 
It  represents  the  arc  which  connects  a  node  with  its  predecessor.   If  a 
node  has  no  predecessor  (i.e.,  a  single  root  node),  then  IVAR  is  set  to 
n  +  1  (the  number  of  arcs  +  1).  During  the  pivot,  IVAR(j)  is  set  to  the 
index  of  the  incoming  arc.  As  the  predecessor  relationship  of  the  j-stem 
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is  reversed,  IVAR  of  each  j-stem  node  must  be  updated  accordingly.  The 
reversal  of  the  j-stem  arcs  corresponds  to  the  reordering  of  the  rows  and 
columns  of  the  current  basis  to  obtain  a  triangular/nearly  triangular 
basis.   IVAR  must  also  exhibit  this  reordering  of  the  basis  elements. 
To  summarize  the  pivot  steps  for  kr(i,  j)  and  k,  (v,  w) : 

(0)  Determine  if  a  new  cycle  is  to  be  formed  (predicted  by  the  ratio 
test:  the  leaving  arc  lies  above  the  join).   If  a  new  cycle  is  formed, 
compute  the  cycle  factor  and  the  dual  of  i,  storing  the  results  in  FAC 
and  U( i)  respectively. 

(1)  Bring  kr(i,  j)  into  the  basis  by  setting  P(j)  <==  i  (with 
appropriate  sign  indicating  arc  orientation);  IVAR(j)  <==  entering  arc, 
and  set  IT  (last  of  the  right  preorder-successors  of  i)  <==  j. 

(2)  Iterate  the  j-stem,  which  includes  the  join.  As  the  j-stem  is 
iterated,  the  predecessor  relationship  is  reversed,  with  P  and  IVAR 
updated  accordingly.   If  the  j-stem  is  on  a  cycle,  the  depths,  D,  of  the 
stem  nodes  are  set  to  zero.  Otherwise,  the  depth  of  j  is  D(i)  +  1  and 
the  preorder-successors  of  j  are  assigned  0(s)  =  D(P(s))  +  1.  The  dual 
of  j  is  computed  as: 

U(j)  =  C(k)  +  MUL(k)  *  U(P(i))        if  P(j)  >  0  or 
U(j)  =  (C(k)  -  U(-  P(i)))/(-  MUL(k))   if  P(j)  <  0; 
where  k  is  the  entering  arc. 

The  duals  of  the  rest  of  the  j-stem  and  their  preorder-successors 
are  computed  similarly. 

If  the  minimum  ratio   is  non-zero,   X   is   also  updated   at   this 
point,   using  equation   (1).      Each  j-stem  node   is  visited    in  turn  by  using 
the  predecessor  relationship   (updating  this  relationship   as    it    is 
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traversed)    and  the  preorder-successors  of  each  j-stem  node  are  visited 
using  the  preorder  relationship,   updating   IT  appropriately. 

(3)     The  i-stem  is  traversed   if  a  new  cycle   is  formed.     During  this 
traversal,  the  dual   variables   are  updated  (replacing  the  reverse  predecessor 
path  temporarily  stored  here)   as  well   as  updating  the  values  of  X.     As  each 
i-stem  node   is  visited,   the  left   and  right  preorder-successors  of  the 
i-stem  node  are  visited   in  preorder,   updating  D,    IT,   and  U.      IVAR  entries 
for  the   i-stem  and   its  preorder-successors  remain  unchanged. 

D.      FURTHER  CONSIDERATIONS 

Priceout,   ratio  test,    and  pivot  provide  the  mechanism  to  move  from 
one  basic  solution  to  a  "better"  basic  solution.     Primal    Simplex  methods 
must  be  designed  to  seek  feasibility  as  well   as  optimal ity.     The  two  most 
common  methods  of  achieving  an  optimal,  basic,   feasible  solution   are  the 
Big-M  method  and  the  two-phase  simplex  method.     The  Big-M  method  uses 
artificial    arcs  with  very  high  costs  to  satisfy  feasibility  initially  and 
the  solution  proceeds  from  this  costly  artificial    start.      In  essence,   the 
Big-M  costs  dominate  the  model    costs  ensuring  that   an  optimal    solution 
will   have  minimal    flow  on   artificial    arcs. 

The  two-phase  method  first   solves  a  related  problem  with   the  same  set 
of  constraints  whose  objective   is   to  minimize  the  sum  of  the  flow  on 
artificial    arcs.      If  an  optimal    solution    is  found  to  the  phase  1  problem 
with   an  objective  function  value  of  zero,   then   all    arcs   pricing  non-zero 
are  eliminated  from  further  consideration   in  phase  2.     The  phase  2 
objective  function,   which    is  the  original   objective  function  of  the 
model,    is   introduced   and  the  optimal    solution   is  sought. 
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An  all  artificial  start  proceeds  by  assigning  IT(i)  =  i,  P(i)  =  i, 
D( i)  =  0  for  i  =  1,  . . .  ,  m  nodes,  and  IVAR  (i)  =  n  +  1  for 
i  =  1,  ...  ,  m,     £ach  node  is  conceptually  assigned  an  artificial  arc 
which  satisfies  the  conservation  of  flow  requirement  at  that  node.  The 
initial  basis  for  m   nodes  is  graphically  depicted  as: 


0(£ 


Figure  20.  All  Artificial  Start 

Each  of  these  artificial  arcs  is  conceptually  assigned  a  multiplier 
of  1.  Thus,  the  dual  associated  with  each  node  is  assigned  a  value  equal 
to  the  cost  of  the  artificial  arc.  The  flow  on  each  of  the  artificial 
arcs  is  set  to  satisfy  conservation  of  flow  requirements. 

Demand  nodes  exhibit  a  negative  external  flow;  to  preserve  non- 
negativity  of  the  right-hand  side,  the  following  adjustments  must  be 
made: 


X( 


<==   -  X(i)  Right-hand   side 

<==  -  P(i)  Arc  orientation 


U(i) 


Dual   values 


The  3ig-'M  nethod   assigns   a  jery  large  (3ig-i*^)   cost  to  eacn  of  the 
artificial    arcs  described   above.     These  costs   are  represented  by  the   ini- 
tial  dual   values   assigned  to  each  node.      If  3ig-M  is  chosen   sufficiently 
large,    an  optimal   minimum  cost   solution  will    not    include  these    /ery 
costly  artificial    arcs  with   non-zero  flow   in  the  basis.      The   algorithm 
will   replace  these  artificial    arcs  with  the  less  costly  "real"   arcs  of 
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the  problem.  The  choice  of  the  Big-M  cost  is  important.  It  must  be 
large  enough  to  drive  the  artificial  arcs  out  of  the  basis,  but  not  so 
large  as  to  cause  numerical  (floating  point)  truncation  errors.  An 
initial  choice  of  twice  the  largest  arc  cost  has  proven  to  be  effective 
in  most  cases. 

An  alternative  to  the  Big-M  method  is  the  two-phase  simplex  approach. 
A  temporary  cost  of  1  is  assigned  to  each  artificial  arc  and  a  temporary 
cost  of  0  is  assigned  to  the  remaining  arcs.  The  algorithm  solves  this 
modified  problem  to  attain  an  initial  feasible  solution.  A  minimum-cost 
feasible  solution  will  have  non-zero  flow  only  on  arcs  with  cost  zero. 
If  an  artificial  remains  in  the  basis  with  non-zero  flow  at  the  end  of 
phase  1,  the  problem  is  infeasible.  Once  phase  1  is  complete,  all  arcs 
with  non-zero  reduced  costs  are  eliminated  from  further  consideration, 
the  correct  costs  are  restored  to  the  arcs  and  the  phase  2  problem  is 
solved,  which  is  the  network  flow  problem  of  interest  restricted  to  admit 
only  feasible  basic  solutions. 

It  is  often  found  that  many  basic  arcs  have  no  successors.  This  is 
especially  true  when  dealing  with  problems  that  have  numerous  sinks.  A 
basis  aggregation  enhancement  to  primal  network  algorithms  has  been 
proposed  by  Bradley,  Brown,  and  Graves  [Ref.  6],  and  has  been  adopted  in 
this  implementation  of  the  GENNET  algorithm.  The  pivotal  update  of  the 
various  arrays  represents  much  of  the  work  performed  by  this  algorithm. 
The  dual  variables,  U,  and  the  depth,  D,  of  leaf  nodes  (i.e.,  nodes  with 
no  successors)  are  uniquely  determined  by  the  knowledge  of  the  leaf 
node's  predecessor  and  the  arc  which  connects  the  leaf  node  to  its 
predecessor.  In  consequence,  these  values  may  be  easily  generated 
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when  required  and  need  not  be  updated  by  every  pivot.  Additionally,  the 
preorder-successor,  IT,  of  this  node  need  no  longer  be  maintained.  While 
it  is  certainly  true  that  any  value  may  be  generated  rather  than  updated, 
the  key  is  to  choose  values  which  may  be  easily  restored. 

Brown  and  McBride  [Ref.  7]  employ  an  array  XM  to  facilitate  the 
"aggregation"  of  nodes.  An  "aggregated  node"  is  a  node  with  no  successors 
whose  depth  and  dual  variable  are  not  explicitly  maintained  and  must  be 
generated  when  required.  The  entries  in  XM  indicate  the  number  of 
successors  of  each  node  which  are  not  currently  explicitly  maintained. 
If  XM(P(r))  =  0,  then  node  r  is  not  an  aggregated  node. 

If  an  aggregated  node  is  encountered  during  priceout,  its  dual  is  gener- 
ated based  on  the  dual  of  its  predecessor  and  the  direction,  multiplier,  and 
cost  of  the  connecting  basic  arc.  Similarly,  the  depth  of  an  aggregated 
node  r  is  generated  as  D(r)  =  0(P(r))  +1.  If  an  aggregated  node,  r,  is 
the  origin  or  terminus  of  the  entering  arc,  it  is  "disaggregated"  by  restor- 
ing its  dual  and  depth,  and  decrementing  XM(P(r)),  the  number  of  aggregated 
successors  of  ?{r) .      IT,  for  the  disaggregated  node  r,  is  updated  by: 
storing  the  preorder-successor  of  ?{r),   TEMP  <==  IT  (P(r));  making  r  the 
preorder-successor  of  P(r),  IT  (P(r))  <==  r;  and  making  the  previous  preorder- 
successor  of  P(r)  the  preorder-successor  of  r,  IT  (r)  <==  TEMP. 

The  outgoing  arc  may  isolate  either  its  head  or   tail  node  with  no 
preorder-successors  and  either  node  may  then  be  aggregated. 

Brown  and  McBride  [Ref.  7]  and  Bradley  et  al .  [Ref.  6]  report  signifi- 
cant computational  savings  using  this  aggregated  node  concept,  especially 
when  dealing  with  real-world  problems  involving  few  sources  and  many  sinks. 
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IV.  MICROCOMPUTER  IMPLEMENTATION 

An  important  emphasis  of  mathematical  programming  has  long  been  on 
the  development  of  computer  codes  which  will  solve  large-scale  problems 
efficiently.  The  network  model,  as  a  specialization  of  the  linear 
program,  was  one  of  the  early  breakthroughs  in  this  area.  As  demonstrated, 
the  unique  character  of  the  network  model  allows  for  efficient  data 
storage  techniques  and  greatly  simplifies  the  computational  requirements 
of  the  simplex  method.  The  obvious  result  is  that  much  larger  problems 
can  be  solved  using  a  network  formulation,  rather  than  a  linear  program, 
on  the  same  computer;  solution  times  and  hence  computational  cost  are 
reduced.  As  computer  memory  becomes  cheaper  and  mainframe  computer 
central  processing  units  (CPU's)  become  faster,  the  problem  size  and  speed 
emphasis  of  mathematical  programming  will  perhaps  diminish.  Indeed, 
problems  of  more  than  one  million  variables  have  already  been  solved 
using  a  model  introduced  by  Geoffrion  and  Graves  [Ref.  20]. 

Advances  in  computer  technology  have  permeated  almost  every  aspect  of 
life  in  the  United  States.  Exhaustive  arithmetic  calculations  can  be 
made  by  anyone  with  a  SIO.OO  calculator.  For  less  than  SIOO.OO,  a 
programmable  calculator/computer  may  be  obtained.  These  radical  develop- 
ments have  been  made  possible  by  advances  in  solid  state  electronics  and 
the  advent  of  the  microprocessor.  A  microprocessor  is  a  collection  of 
many  thousand  electronic  logical  gates  created  as  microscopic  circuits  on 
small  pieces  of  silicon,  known  as  chips  [Ref.  21].  Microprocessors  are 
perhaps  best  known  as  the  controlling  device  in  home  and  arcade  video 
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games.  They  also  serve  more  practical  purposes  as  industrial  process 
controllers  and  as  the  CPU  of  the  microcomputer. 

The  most  popular  microprocessors  used  in  today's  small  computers  are 
the  8080  manufactured  by  Intel,  the  Z80  made  by  Zilog,  the  6800  manufac- 
tured by  Motorola,  and  the  6502  developed  by  MOS  Technology  [Ref.  22]. 
Each  of  these  are  8-bit  microprocessors,  indicating  that  the  basic  memory 
unit  is  8  bits  wide  (a  byte).  These  microprocessors  have  16  address 
lines.  Each  of  these  lines  may  be  in  one  of  two  states:  high  or 
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low  (on  or  off;  0  or  1).   As  such,  these  microprocessors  may  address  2  , 
or  65,536,  separate  memory  locations.  A  microcomputer  with  65,536 
addressable  memory  locations  is  known,  somewhat  inaccurately,  as  a  64K  (K 
stands  for  kilobyte)  microprocessor. 

The  addressable  memory  in  a  microprocessor  is  not  all  available  to 
the  user.  Most  microcomputer  systems  reserve  some  of  that  memory  for  the 
use  of  the  monitor,  various  languages,  and  operating  systems.  The  Apple 
II  Plus  microcomputer,  in  its  64K  configuration  running  the  Apple  UCSD 
(University  of  California  at  San  Diego)  Pascal  Operating  System,  has  a 
maximum  space  of  39,900  bytes  for  a  user  program  and  variables  [Ref.  23]. 

In  the  last  two  years,  16-  and  32-bit  microprocessors  have  emerged, 
as  well  as  some  microprocessors  with  20  or  more  address  lines  (making 
them  capable  of  addressing  more  than  one  million  memory  locations) 
[Ref.  21].  Certainly,  size  distinctions  between  microcomputers  and 
minicomputers  have  diminished  and  microcomputers  may  soon  be  challenging 
mainframes  in  many  applications. 

E.  M.  L.  Beale,  in  his  1980  Blackett  Memorial  Lecture  on  the  relation- 
ship between  operations  research  and  computers  [Ref.  24],  recognizes  the 

65 


emergence  of  microcomputers  and  their  potential  as  a  powerful  tool. 
There  has  been  a  plethora  of  software  developed  for  microcomputers  for 
use  in  the  statistical  and  trend  analysis  areas  of  operations  research. 
There  is,  however,  a  dearth  of  mathematical  programming  software  available 
to  the  microcomputer  user  and  surprisingly  little  published  research  in 
this  area. 

The  potential  of  small  computers  has  not  escaped  all  mathematical 
programmers  and  researchers.  In  1979,  a  feasibility  study  [Ref.  25]  was 
conducted  to  explore  the  possibility  of  implementing  mathematical  program- 
ming algorithms  on  minicomputers.  The  study  implemented  two  shortest 
path  algorithms  on  two  different  minicomputers.  It  concluded  that 
minicomputers,  although  slower  than  mainframes,  were  acceptable  vehicles 
for  shortest  path  algorithms.  The  study  also  hypothesized  that  modern 
minimum  cost  flow  network  algorithms  may  also  be  excellent  candidates  for 
minicomputer  implementation. 

Also  in  1979,  F.  P.  Wyman  [Ref.  26]  reported  the  implementation  of  an 
out-of-kilter  algorithm  for  the  solution  of  pure  minimum  cost  network 
flow  models.  Additional  implementations  of  "textbook"  style  PERT,  CPM, 
and  SIMPLEX  codes  have  been  reported  in  hobby  computer  journals  [e.g., 
Ref.  27]. 

R.  H.  Duff  [Refs.  28,  29]  reported  the  development  of  a  comprehensive 
microcomputer-based  network  optimization  package.  Duff's  microcomputer 
package  is  capable  of  solving  pure  minimum  cost  network  flow  problems, 
elastic  network  problems  (in  which  flow  conservation  may  be  violated  for 
a  "price"),  and  nonlinear  network  problems.  The  algorithms  chosen  were 
state-of-the-art  optimization  algorithms  designed  to  minimize  storage 

66 


requirements  and  execution  time.  Algorithms  with  these  characteristics 
are  obvious  candidates  for  implementation  on  microcomputers  where  memory 
is  a  limited  commodity  and  CPU  processing  is  in  the  one  to  three  megahertz 
(MHz)  range. 

To  complete  Duff's  network  optimization  package  and  make  it  the 
most  versatile  network  optimization  package  available  on  anything  but  a 
mainframe  computer  (and  perhaps  not  even  there),  all  that  is  needed  is  an 
efficient  generalized  network  algorithm. 

The  generalized  network  algorithms  presented  by  Jensen  and  Barnes 
[Ref.  1]  and  Kennington  and  Helgason  [Ref  2]  are  thought  to  be  too 
cumbersome  and  inefficient  for  the  "lean"  world  of  microcomputing.  Elam 
et  al .  [Ref.  30]  report  the  development  of  a  fast  and  efficient  general- 
ized network  code,  but  descriptions  of  that  code  are  spread  throughout 
the  literature  [Refs.  5,  30,  31,  32]  and  provide  no  clear  explanation  of 
the  basis  update  mechanism.  As  such.  Brown  and  McBride's  [Ref.  7]  GENNET 
Algorithm,  as  described  in  detail  in  the  previous  sections,  has  been 
chosen  as  the  most  suitable  addition  to  the  microcomputer  network  optimiza- 
tion package.  For  the  sake  of  brevity,  the  microcomputer-based  network 
optimization  package  herein  described  will  be  referred  to  as  Micronet. 

The  host  computer  for  Micronet  is  an  Apple  II  Plus  with  54K  of 
memory.  This  is  certainly  not  the  most  powerful  microcomputer  (with  some 
competitors  featuring  addressing  capabilities  of  15  megabytes  and  running 
at  speeds  of  8  MHz),  but  it  is  one  of  the  most  popular  with  a  population 
of  more  than  400,000  units  [Ref.  33 J. 

It  is  important  to  recognize  that  the  host  machine  is  not  at  issue 
here.  This  project  was  begun  in  1979  and  the  Apple  II  was  representative 
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of  the  technology  emerging  at  that  time.  The  Apple  II  has  continued  to 
be  an  extremely  popular  and  representative  8-bit  microcomputer.  The 
purpose  of  Micronet  is  to  explore  the  capabilities  and  suitability  of  a 
microcomputer  as  a  tool  of  mathematical  programming.  Software  design 
techniques  used  in  this  optimization  package  are  applicable  to  other 
microcomputer  systems  as  well.  Larger  and  faster  microcomputers  will  only 
enhance  the  capabilities  of  software  systems  such  as  Micronet,  making 
mathematical  programming  on  a  microcomputer  not  only  a  viable,  but  an 
attractive  option. 

In  this  spirit,  the  Apple  II  Plus  continues  to  be  used  as  the  host 
computer  for  the  further  development  of  microcomputer-based  network 
codes.  The  programming  language  used  is  UCSD  PASCAL.  This  language  is 
fast  becoming  one  of  the  most  popular  microcomputer-based  high  level 
computer  languages.  Certainly  BASIC  ranks  as  the  most  popular,  however, 
versions  of  BASIC,  both  interpreted  and  compiled,  lack  standardization. 
UCSD  PASCAL,  while  all  implementations  are  not  identical,  provides  a  more 
standardized  vehicle  for  the  development  of  microcomputer  programs.  Duff 
[Ref.  28]  discusses  the  choice  of  PASCAL  in  great  detail;  his  reasoning 
remains  sound  and  will  not  be  repeated. 

It  is  appropriate,  however,  to  discuss  some  limitations  of  the  host 
microcomputer  and  the  Apple  implementation  of  UCSD  PASCAL;  these  limita- 
tions profoundly  impact  the  design  of  any  microcomputer  mathematical 
programming  code. 

Pure  network  problems  are  characterized  by  columns  of  the  constraint 
matrix  which  have  non-zero  entries  of  only  +1  or  -1.  This  unique 
characteristic  eliminates  the  need  for  floating  point  arithmetic  and 
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eliminates  the  requirements  for  multiplication  and  division.  As  such, 
mathematical  precision  is  not  an  issue  in  pure  network  codes.  Generalized 
networks  possess  a  similar  structure  with  each  column  of  the  constraint 
matrix  having  a  +1  and  a  multiplier  -m.  ,  where  -m.  may  be  any 
floating  point  number.  The  optimal  flows  in  a  generalized  network  are 
therefore  not  necessarily  integer,  thus  floating  point,  as  opposed  to 
integer,  arithmetic  is  required. 

The  Apple  II  Plus  implementation  of  UCSD  PASCAL  provides  for  a  32-bit 
representation  of  floating  point  numbers  [Ref.  34],  A  representation 
using  24  bits  for  the  mantissa  of  a  floating  point  number  and  8  bits  for 
the  exponent  provides  six  or  seven  significant  figures  of  precision  with 
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a  dynamic  range  of  10"   to  10  ^     [Ref.  35].  This  precision  is 
roughly  equivalent  to  IBM  360/370  single  precision.  Although  this  is  a 
limitation  of  Apple  PASCAL,  it  is  not  considered  to  be  severely  debilitating 
as  applied  to  generalized  networks.  The  purpose  of  the  multiplier 
(-m,  )  is  to  model  the  transformation  of  units  of  flow  or  change  the 
commodity  amount  as  it  flows  through  an  arc  [e.g.,  Ref.  5].  These  purposes 
of  arc  multipliers  can  usually  be  accommodated  with  two  or  three  significant 
figures.  As  such,  an  arithmetic  precision  level  of  10"  has  been 
chosen  for  use  in  this  implementation  of  GENNET. 

UCSD  PASCAL  requires  static  dimensioning  of  arrays.   Proqramminq 
techniques  which  dynamically  allocate  memory  at  execution  time  based 
on  problem  size  and  type  cannot  be  used  to  provide  better  memory  management. 
This  is  not  a  critical  economic  issue  on  a  dedicated  microcomputer,  as  it 
is  on  a  mainframe  computer,  but  it  does  unnecessarily  limit  problem  size 
and  flexibil ity. 
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Source  statements  are  compiled  into  a  standard  UCSD  pseudo-code 
(p-code)  which  is  then  interpreted  at  execution  time.  This  system  allows 
for  a  standard  p-code  across  all  machines  and  only  requires  a  machine- 
dependent  run  time  interpreter.  Execution  speed  is  therefore  not  as  fast 
as  would  be  expected  from  a  truly  compiled  language,  but  not  as  slow  as 
an  interpreted  BASIC.  The  most  disturbing  feature  of  Apple  II  PASCAL, 
when  used  as  a  vehicle  for  mathematical  programming,  has  proven  to  be 
compilation  speed.  The  nominal  compilation  rate  for  Apple  II  UCSD  PASCAL 
is  two  hundred  source  statements  per  minute.  This  means  that  a  large 
program  might  take  20  to  30  minutes  to  compile  and  that  requirement  can 
be  quite  annoying  when  a  program  is  in  the  developmental  stage. 

Other  less  significant  limitations,  which  are  included  for  the  sake 
of  completeness,  are  the  editable  file  size  and  maximum  procedure  size. 
The  operating  system  editor  accommodates  files  of  sizes  up  to  40  blocks 
[Ref.  23].  A  block  represents  500  bytes  of  information  stored  on  an 
Apple  II  floppy  diskette.  To  compose  and  edit  large  programs,  several 
text  files  must  be  created  and  connected  together.  The  maximum  size  of  a 
compiled  procedure  or  function  is  1,200  bytes  (plus  any  local  variables). 
This  requires  programs  to  be  broken  into  smaller  "pieces"  than  may  be 
logically  convenient.  These  latter  limitations  are  obviously  not  major 
faults,  but  are  interesting  "quirks,"  which  must  be  contended  with. 

Once  the  Apple  PASCAL  operating  system  is  resident  in  memory,  less 
than  40K  bytes  are  available  for  user  programs  and  variables  [Ref.  23]. 
As  such,  any  mathematical  programmer  must  be  extremely  conscious  of  the 
classical  space-speed  tradeoffs.  Certainly  the  efficient  data  storage 
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techniques  and  basis  update  mechanisms  exhibited  by  the  GENNET  [Ref.  7] 
algorithm  are  mandated. 

The  UCSD  PASCAL  system  allows  "segmentation"  or  "overlaying"  of 
program  components  for  large  programs.  When  a  portion  of  code  is  no 
longer  required,  it  may  be  "swapped"  out  of  memory  and  replaced  by  another 
segment  of  code.  To  maximize  the  amount  of  memory  available  for  problem 
representation,  it  is  desirable  to  have  as  little  memory  occupied  by 
program  code  as  possible.  However,  as  more  program  overlaying  is  done, 
more  disk  accesses  are  required  by  code  swapping.  Disk  accesses  are  slow 
and  the  programmer  is  again  presented  with  the  ubiquitous  space-speed 
tradeoff. 

Micronet  consists  of  three  major  components:  a  master  program  or 
driver,  an  editor,  and  a  solution  module.  These  components  all  have 
access  to  a  system  library  containing  often-used  functions.  An  online 
use  manual,  or  "help"  feature,  is  eventually  envisioned  for  Micronet  but 
has  not  been  completely  implemented  as  yet.  The  solution  module  has 
several  submodules  for  solving  the  various  types  of  network  flow  problems: 
pure  minimum  cost  (GNET),  nonlinear  (NLPNET),  elastic  with  fixed  charges 
(ENET),  and  generalized  networks  (GENNET).  A  conceptual  view  of  .Micronet 
is  presented  in  Figure  21. 

Duff  [Ref.  28]  gives  a  detailed  description  of  organization,  use, 
and  the  characteristics  of  each  component  of  Micronet  (with  the  exception 
of  GENNET),  The  driver,  editor,  and  solution  module  have  been  appropri- 
ately modified/amended  to  allow  Micronet  to  solve  generalized  networks. 
A  brief  description  of  Micronet  will  be  given  here   for  continuity. 
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MICRONET  ORGANIZATION 


EDITOR 
MODULE 


GNET 


DRIVER 


SOLUTION 
MODULE 


->    NLPNET 


ENET 


ON  LINE 

USE 
MANUAL 


GENNET 


Figure  21.  ^icronet  Organization 

The  iiiaster  program  coordinates  the  operation  of  the  optimization 

package  oy  passing  control  to  the  editor  for  data  nanioulation,  to  the 

solution  ncdule  to  solve  a  network  model,  or  to  the  ^ASCAL  ooerating 

system  to  exit  .Micronet.  Upon  completion  of  work  by  the  solution  or 

editor  module,  control  is  passed  back  to  the  driver  program.   The  driver 

is  chained  to  the  solution  module  and  editor  so  that  when  control  is 

passed  between  one  of  the  three  main  package  comoonents,  only  that 
component  is  in  memory. 
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The  purpose  of  the  editor  is  to  create,  alter,  transfer,  and  browse 
data  files.  Data  files  are  essentially  network  problem  description  files 
consisting  of  a  header  record  and  records  for  nodes  and  arcs.  Nodes  may 
be  specifically  assigned  attributes  (e.g.,  flow  bounds,  external  flow 
requirements),  or  the  problem  may  be  described  using  artificial  arcs 
representing  the  external  flow  characteristics  of  sources  and  sinks. 
Problem  files  may  be  created  from  the  keyboard  or  read  from  a  text  file 
in  SHARE  format  [e.g.,  Ref.  36].  Examples  of  the  three  types  of  records 
follow: 


Header  Record 


Problem  Name: 
Problem  Type: 
Number  of  Nodes 
Number  of  Arcs: 
Date  Created: 


GENTRANS 

GENERALIZED 

15 

30 

1  Oct  81 


Date  Last  Updated:  15  Nov  81 


Arc  Record  (A  Generalized  Arc) 


Arc  Name: 
Source  Node: 


Destination  Node:  5 


Lowerbound: 

2.0 

Upperbound: 

100.0 

Initial  Flow: 

0 

Unit  Cost: 

10.50 

Multipl ier : 

1.18 
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Node  Record 
Node  Name:        Phoenix 
Node  Number:      5 
Node  Kind:        Demand 
Net  Flow:         3.0 
Lower  Range:      0 
Upper  Range:      0 
Lower  Penalty:     0 
Upper  Penalty:     0 
The  ranges  and  penalties  indicated  on  the  node  record  are  applicable  to 
elastic  programming  and  are  discussed  by  Duff  [Ref.  28]. 

As  the  information  is  entered  into  the  problem  record,  it  is  screened 
for  consistency.  Once  created,  a  file  can  be  altered,  browsed,  removed, 
or  transferred  by  choosing  the  appropriate  editor  option. 

If  the  solution  module  is  chosen  from  the  command  level,  the  user  is 
prompted  to  insert  a  disk,  containing  the  problem  file,  into  one  of  the 
Apple  II's  disk  drives.  Micronet  then  displays  a  catalog  of  the  disk  and 
the  user  selects  the  problem  to  be  solved.  At  this  point,  the  header 
record  is  read  and  based  on  the  problem  type,  the  appropriate  solution 
submodule  is  automatically  invoked.  Figure  22  displays  the  solution  path 
selection  logic. 

If  the  problem  type  is  a  generalized  network,  the  GENNET  solution 
module  is  chosen.  This  module  implements  the  previously  described  GENNET 
Algorithm  in  a  segmented  PASCAL  Program.  The  program  is  split  into  four 
segments  as  shown  below: 

Main  Program 

Input  -^^    Initialization  ^  Solution  ^  Report 
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The  segmentation  of  the  program  represents  this  space-speed  trade-off 
decision  that  must  be  made.  The  main  program  calls  each  of  the  segments 
in  turn,  "swapping"  out  of  memory  the  previously  used  segment.  The  main 
program  declares  all  global  variables  and  thereby  provides  for  communica- 
tion between  segments.  These  variables  remain  in  memory  throughout  the 
use  of  the  solution  module.  Variables  local  to  the  individual  segments 
occupy  memory  only  when  that  segment  is  active. 

The  input  segment  first  presents  the  user  with  a  menu  requesting 
a  destination  for  the  results  (disk  file,  printer,  serial  interface,  or 
terminal),  whether  a  pivot-by-pivot  trace  is  required,  and  whether  a 
listing  of  the  problem  arcs  and  nodes  is  desired.  The  input  segment 
establishes  the  maximum  problem  size  at  100  nodes  and  500  arcs.  The 
input  file  is  read,  and  problem  variables  are   initialized.  Lower  bounds 
are  translated  out  and  the  resulting  cost  of  the  lower  bound  modification 
is  recorded. 

During  this  phase  of  the  solution  procedure,  the  user  is  required  to 
interact  with  the  computer  by  answering  menu-driven  questions,  inserting 
a  problem  disk,  and  choosing  a  problem.  Great  care  has  been  exercised  to 
ensure  that  an  incorrect  user  response  will  illicit  a  helpful  (or  admonish- 
ing) message  from  the  program  rather  than  a  premature  program  termination. 
Having  read  the  problem,  the  header  record,  arc  records,  and  node  records 
have  been  converted  to  the  streamlined  data  structures  of  GENNET. 

During  the  initialization  segment,  a  Big-M  start  is  set  up,  with  the 
initial  basis  appearing  as  a  forest  of  one-trees.  Big-M  is  determined  by 
doubling  the  largest  arc  cost  of  the  problem.  The  arcs  are  sorted  by 
head  node  and  stored  in  that  order  in  the  array  T.  T(k)  contains  the 
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c 


LINEARIZE 


GMET 
START 


♦   USER  iNTERACTlON  REQUIRED 
♦*    MAY   BE  ELASTIC  OR   "O-U" 
+    WARNING  MESSAGE  ISSUED 


Figure  22.  Solution  Module  Selection  Logic 
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tail  of  arc  k.  The  head  node  pointers  in  H  are  also  set  up  at  this  time 
and  the  arrays  describing  arc  characteristics  are  sorted  corresponding  to 
the  order  of  T.  Non-negativity  of  the  right-hand  side  is  enforced  at  this 
point,  with  the  dual  variables  and  predecessors  being  appropriately  updated 

Having  initialized  the  problem,  the  program  now  executes  the  solution 
segment.  Two  new  arrays  must  be  introduced  at  this  point  for  operation 
of  the  candidate  queue,  making  this  the  most  space-critical  segment.  As 
such,  the  coding  in  this  segment  is  terse  to  maximize  efficiency  and 
acceptable  problem  size.  Larger  problems  could  be  accommodated  if  the 
solution  module  were  segmented  into  smaller  components  such  as  priceout, 
ratio  test,  and  pivot.  This  would,  however,  cause  multiple  memory-to-disk 
"swaps"  for  each  pivot,  and  disk  operations  are  extremely  slow  compared 
with  core-memory  operations.  The  decision  has  been  made  to  sacrifice 
problem  size  in  favor  of  solution  speed.  To  perform  memory  "swaps"  for 
each  pivot  is  considered  to  be  prohibitively  slow. 

The  operation  of  the  solution  module  is  split  into  the  three  logical 
steps  of  the  simplex  method:  priceout,  ratio  test,  and  pivot.  The 
candidate  queue  of  "interesting  nodes  and  good  arcs"  is  initialized  with 
the  sink  nodes--if  none  are  found,  the  queue  may  be  initialized  with  any 
node.  Completion  of  the  solution  module  is  determined  by  exceeding  a 
preset  maximum  number  of  pivots  or  by  pricing  out  every  non-basic  arc  and 
determining  that  none  price  favorably. 

If  the  pivot  count  has  not  been  exceeded  and  the  solution  module  ends 
with  a  possible  optimal  solution,  the  final  report  segment  is  called. 
This  final  segment  first  determines  if  there  are  any  artificial  arcs  left 
in  the  basis  with  non-zero  flow.   If  there  are,  the  solution  is  not 
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feasible.  The  sum  of  the  flow  on  the  artificials  is  recorded  and  the 
problem  is  resolved  using  a  larger  Big-M  cost.  After  this  solution 
iteration,  if  a  non-feasible  solution  is  again  obtained,  the  total  flow 
on  the  artificial  arcs  is  compared  to  the  previous  total.   If  the  arti- 
ficial flow  has  decreased,  Big-M  is  again  increased  and  the  problem  is 
again  resolved.   If  the  artificial  flow  has  not  decreased  from  the 
previous  iteration,  the  solution  process  is  terminated  and  the  problem  is 
declared  infeasible. 

Figure  23  displays  the  input  arrays  of  a  small  generalized  network 
problem  and  Figure  24  exhibits  the  report  of  an  optimal  solution.  The 
final  flows  indicated  on  the  arcs  are  flows  in  addition  to  arc  lower 
bounds. 
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APPLE^4ET  -  GENNET  MODULE 
I  OF  MAR  3  2 


(500  ARC  130  NODE  VERSION] 
DATE:  19  APR  3  2 


FILE:  PR0B:GT2.NET 
NUMBER  OF  NODES  =  15 


CREATED:  13  SEP  81 
NUMBER  OF  ARCS  =  30 


UPDATED:  14  APR  82 


ARC  LIST  . 

•  « 

ARC 

FROM 

TO 

UNIT 

UPPER 

LOWER 

NAME 

NODE 

NODE 

COST 

GAIN 

BOUND 

BOUND 

1 

4 

3 

33.84 

0.99 

1003.00 

0  .00 

2 

2 

3 

15.47 

1.00 

1030.33 

0.30 

3 

1 

5 

53.54 

0.74 

1000.00 

3.33 

4 

2 

5 

26.75 

0.74 

1000.03 

3.00 

5 

3 

5 

73.49 

1.00 

1003.00 

0.00 

6 

5 

5 

52.52 

1.00 

1303.03 

0.30 

7 

3 

6 

35.12 

0.91 

1030.00 

0.33 

8 

5 

6 

11.12 

1.30 

1000.00 

0.00 

9 

4 

7 

59.55 

1.17 

1300.00 

3.30 

10 

2 

7 

88.38 

1.35 

1000.30 

2.30 

11 

4 

8 

84.12 

1.00 

1030.03 

3.03 

12 

2 

8 

21.85 

0.92 

1300.30 

3.33 

13 

4 

9 

3.46 

1.00 

1303.30 

0.30 

14 

3 

9 

29.72 

1.00 

1330.03 

0.30 

15 

4 

10 

6.12 

1.00 

1333.00 

3.03 

16 

2 

10 

31.08 

0.96 

1330.03 

3.00 

17 

3 

10 

1.07 

1.07 

1003.00 

0.03 

18 

5 

10 

44.44 

1.30 

1300.03 

0.30 

19 

1 

11 

67.15 

0.91 

1003.00 

0.00 

20 

2 

11 

59.33 

0.79 

1030.30 

3.03 

21 

3 

11 

50.45 

1.17 

1033.30 

0.00 

22 

5 

11 

71.42 

1.30 

1303.30 

0.00 

23 

2 

12 

8.88 

1.18 

1003.00 

0.00 

24 

1 

13 

28.22 

0.33 

1003.00 

0.33 

25 

4 

13 

77.34 

1.00 

1000.30 

0.00 

26 

3 

13 

45.50 

1.00 

1333.30 

3.00 

27 

5 

13 

20.57 

0.83 

1003.00 

3,00 

28 

4 

14 

37.75 

1.13 

1333.03 

3.30 

29 

2 

14 

18.16 

0.98 

1300.00 

0.30 

30 

3 

15 

67.52 

1.00 

1330.30 

0.33 

NODE  LIST 

...  [  FOR 

THOSE  : 

:J0DES  EXPLICIT! 

LY  IN  THE  DATA 

FILE  ] 

NODE 

NODE 

NET 

NODE 

NAME 

NUMBER 

FLOW 

STATUS 

1 

1 

22, 

,86 

SUPPLY 

2 

2 

177, 

,14 

SUPPLY 

6 

5 

-19. 

.39 

DEMAND 

7 

7 

-3. 

,64 

DEMAND 

8 

8 

-24, 

,92 

DEMAND 

9 

9 

-9, 

,33 

DEMAND 

10 

13 

-14. 

,07 

DEMAND 

11 

11 

-56. 

,91 

DEMAND 

12 

12 

-2. 

.45 

DEMAND 

13 

13 

-30, 

,93 

DEMAND 

14 

14 

-21. 

,75 

DEMAND 

15 

15 

-16. 

,55 

DEMAND 

Fi 

gure  23 

.  General ized 

Network  Problem 
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OPTIMALITY  OBTAINED  IN  23  PIVOTS 


ARC 

FROM 

TO 

COST 

UP-50UND 

DUAL 

FLOW  ABOVE 

24 

1 

13 

23.2203 

1000.00 

-49.7936 

22.3500 

2 

2 

3 

15.47Cf! 

1000.00 

-32.9223 

115.039 

7 

3 

6 

35.1200 

1000.00 

-48.3923 

15.4126 

ART 

4 

4 

O.BnD0O 

INF 

-220.000 

0.00000 

4 

2 

5 

26.7500 

1000.00 

-80.5517 

7.24931 

3 

5 

6 

11.1200 

1000.00 

-91.7718 

5.35449 

10 

2 

7 

33.3S00 

1000.30 

-114.435 

1.4  3  395 

12 

2 

8 

2i.an00 

1003.00 

-59.5450 

27.0870 

14 

3 

9 

29.7200 

1000.33 

-73.1123 

9.33303 

17 

3 

10 

1.07  0:^0 

1000.00 

-46.2254 

13.1495 

21 

3 

11 

50.4530 

1300.00 

-34.4891 

48.5413 

23 

2 

12 

8.88000 

1000.03 

-35.4257 

2.37527 

26 

3 

13 

45.6000 

1000.00 

-93.9923 

11.9562 

29 

2 

14 

18.1600 

1000.03 

-52.1243 

22.2041 

30 

3 

15 

57.6230 

1003.30 

-115.012 

15.5530 

VALUE  OF  THE  OBJECTIVE  FUNCTION=  3.94934E3 


FINAL  DUMP 


PIVOTS: 23   DEGENERATE  PIVOTS: 3   CYCLES  CREATED: 4 


Figure  24.  Generalized  Network  Solution 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 

This  project  was  initiated  for  two  reasons.  The  first  was  to  better 
understand  generalized  networks  and  to  gain  an  appreciation  for  and 
familiarity  with  the  data  structures  and  design  of  large-scale  mathemati- 
cal programming  algorithms.  Secondly,  an  objective  was  to  explore  the 
suitability  of  a  microcomputer  as  a  tool  of  operations  research. 

The  design  and  implementation  of  a  large-scale  mathematical  program- 
ming project  has  been  presented  here  in  great  detail  (only  the  memory 
size  of  the  Apple  II  microcomputer  limits  the  code  to  medium-sized 
problems).  In  mathematical  programming  in  general  and  microcomputer 
programming  in  particular,  the  requirement  to  use  sparse  data  structures 
and  efficient  computational  mechanisms  cannot  be  overemphasized.  The 
programmer  must  vigorously  search  for  ways  to  condense  coding  segments 
and  use  mathematical  simplification/insight  to  reduce  object  code  and 
computational  overhead.  At  the  same  time,  a  computer  program  must  have 
an  easily  accommodated  user  Interface  if  it  is  to  be  of  real  value. 

While  the  programmer  must  be  terse  and  exceedingly  efficient  when 
designing  a  solution  module,  that  same  programmer  must  be  lavish  when 
designing  the  man-machine  interface.  The  primary  communications  interface 
for  the  microcomputer  is  the  keyboard.  Mis-stroked  keys  must  be  screened, 
and  incorrect  user  inputs  must  be  tolerated  by  the  program.  Single 
stroke  responses  to  menus  have  been  found  to  be  the  best  method  of 
communication.  When  numerical  data  is  required,  the  user  must  be  provided 
the  luxury  of  easily  amending  the  input. 
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The  network  optimization  package,  as  described  in  this  thesis, 
accommodates  both  the  lean  solution  module  and  forgiving  user  interface 
requirements.  The  editor  component  of  the  package  provides  the  primary 
man-machine  interface  and  is  designed  to  prompt  and  edit  all  inputs.  The 
solution  modules  have  been  coded  very  succinctly  and  utilize  the  extremely 
efficient  data  structures  of  which  the  GENNET  algorithm  is  exemplary. 

This  project  has  demonstrated  that  serious  mathematical  programming 
can  be  accommodated  by  microcomputers.  The  primary  drawback  of  micro- 
computers in  this  endeavor  is  considered  to  be  the  slow  compilation 
rate.  However,  the  disadvantages  associated  with  slow  compilation 
represent  an  initial  development  cost  and  as  such  are  considered  acceptable, 
given  the  utility  of  the  resulting  product. 

The  integrated  structure  of  Micronet  eases  the  burden  of  data  input 
requirements.  Laborious  keyboard  sessions  may  be  avoided  if  problems, 
existing  as  textfiles  on  mainframes,  are  transferred  to  a  microcomputer 
textfile.  The  Micronet  editor  has  the  ability  to  create  a  problem  file 
of  arc  and  node  records  from  a  textfile  in  SHARE  format  [e.g.,  Ref.  36]. 
Alternatively,  a  computer-based  problem  generation  technique  could  be 
employed.  In  so  doing,  a  microcomputer  (or  any  computer)  is  capable  of 
solving  problems  which  are  much  larger  than  can  be  reasonably  created 
from  a  keyboard.  Real -life  problems  to  be  solved  by  microcomputers  may 
well  be  based  on  microcomputer  data  files.  In  this  case,  data  input 
requirements  could  again  be  automated. 

Microcomputers  continue  to  evolve  at  an  extremely  rapid  rate. 
Prices  are  falling  while  memory  addressing  capability  and  CPU  speed  are 
being  enhanced.  Some  microcomputers  can  address  up  to  16  megabytes  of 
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memory,  perform  double  precision  arithmetic  using  a  64-bit  representation 
of  floating  point  numbers,  and  have  internal  clock  speeds  eight  times  as 
fast  as  the  Apple  II.  As  such,  it  is  considered  pointless  to  discuss 
solution  speeds,  precision,  and  problem  size.  The  current  code  accommo- 
dates problems  of  sizes  up  to  100  nodes  and  500  arcs,  with  execution 
times  averaging  .5  seconds  per  simplex  iteration.  Certainly  any  published 
figures  could  easily  be  eclipsed  by  new  entries  in  the  microcomputer 
market,  costing  little  more  than  the  original  list  price  of  a  fully 
equipped  Apple  II.  This  project  has  demonstrated  that  efficient  algorithms 
can  be  constructed  and  implemented  on  microcomputers  for  the  broad  class 
of  problems  which  may  be  modelled  by  a  generalized  network.  The  sparse 
data  structures  and  computational  efficiency  afforded  by  modern  mathematical 
programming  codes  are  well-suited  for  implementation  on  microcomputers. 

Managers,  scientists,  small  businesses,  and  government  agencies 
purchasing  microcomputers  for  other  than  mathematical  programming  purposes 
can  feasibly  add  a  microcomputer  optimization  package.  Small  colleges 
and  businesses  often  cannot  afford  access  to  a  mainframe  computer  optimi- 
zation package.  The  development  of   microcomputer-based  mathematical 
programs  make  these  management  and  decision-making  tools  available  to  a 
much  wider  audience. 

Mathematical  programming  software  for  microcomputers  is  currently  in 
scarce  supply,  while  the  population  and  availability  of  microcomputers  is 
growing  at  a  breathtaking  rate.   As  managers  and  operations  researchers, 
we  must  take  full  advantage  of  the  power  and  flexibility  of  the  micro- 
computer and  begin  to  export  mathematical  programming  methods  to  the 
large  audience  of  microcomputer  users.   It  is  doubtful  that  mainframe 
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computers  will  be  replaced  as  the  primary  tool  of  the  operations  researcher 
and  mathematical  programmer.  However,  the  advantages  and  techniques  of 
scientific  management  through  mathematical  programming  should  be  made 
available  to  the  "common  man"  (or  at  least  the  common  manager)  through 
the  vehicle  of  microcomputers. 
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